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ABSTRACT 


Power  electronics  have  advanced  to  the  point  that  they  can  be  considered  for  use 
in  large,  high  power  dc  electrical  distribution  systems.  The  proposed  Navy  DC  Zonal 
Electrical  Distribution  System  (DC  ZEDS)  being  designed  for  the  new  DD-21  utilizes 
a  rectified  ac  generator  ou^ut  which  is  filtered  and  stepped  to  usable  voltages  by 
local  dc-dc  converters.  One  characteristic  of  the  high-bandwidth  converters  is  a 
negative  input  impedance,  which  when  coupled  wifii  an  LC  input  filter,  can  lead  to 
system  instabilities.  This  thesis  examines  various  stability  criteria  for  determining 
parameters  of  the  dc  link  filter.  Comparisons  between  a  simplified  system  model,  a 
model  using  subsystem  impedances  and  an  Advanced  Continuous  Simulation 
Language  (ACSL)  model  of  a  reduced-order  system  are  made.  Simulations  were 
conducted  to  verify  the  vaUdity  of  the  stability  criteria.  The  ACSL  model  provides 
an  extremely  useful  tool  in  evaluating  the  response  of  various  system  parameters  to 
changes  in  design  values.  The  design  criteria  examined  in  this  thesis  can  be 
ultimately  used  to  provide  design  specifications  to  future  vendors. 
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I. 


PROJECT  BACKGROUND 


Many  recent  issues  have  evolved  which  warrant  the  need  and  provide  the 
motivation  for  a  new  way  of  thinking  about  electric  power  distribution  on-board  Navy 
ships.  In  recent  years,  the  Navy  has  pushed  for  more  use  of  Commercial-Off-The-Shelf 
(COTS)  parts  in  their  new  ship  designs  and  existing  ship  upgrades.  Hand-in-hand  with 
this  trend  is  the  updating  and  designing  of  systems  using  the  latest  and  most  advanced 
technology  available  on  the  world  market.  Budgets  are  currently  limited,  which  is 
forcing  the  need  for  cost  effective  designs.  As  ships  are  getting  more  complex,  the 
possibility  of  automated  systems  is  growing.  The  move  towards  automation  is  also 
driven  by  the  economics  of  streamlining  the  manning  of  future  Navy  ships.  Naval 
architecture  is  experimenting  with  unique  designs  which  require  iimovative  ship  system 
designs  to  allow  for  more  flexibility  in  their  arrangements.  Currently,  the  Navy 
implements  a  radial  distribution  configuration  for  electric  power  on-board  its  vessels.  A 
suggested  alternate  option  is  that  of  zonal  distribution. 

The  Direct  Current  Zonal  Electrical  Distribution  System  (DC  ZEDS)  is  a 
component  of  the  Integrated  Power  System  (IPS)  currently  being  proposed  for  use  on  the 
new  DD-21  design.  The  basic  concept  for  the  IPS  is  to  use  a  common  set  of  prime 
movers  to  generate  power  for  both  propulsion  and  distribution.  It  is  fundamentally 
different  from  the  radial  design  in  many  respects.  These  differences  provide  numerous 
advantages  as  outlined  below. 

A.  RADIAL  DISTRIBUTION 

The  radial  electric  power  distribution  concept  involves  distributing  60  Hz 
alternating  current  (ac)  power  produced  by  a  limited  number  of  generators  to  loads 
throughout  the  ship  in  a  pyramid  type  configuration.  Branching  occurs  from 
switchboards  which  route  power  to  power  panels  located  in  the  general  area  of  the  loads 
that  they  supply.  Loads  which  are  considered  non-vital  receive  power  from  a  single 
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source.  Loads  which  are  considered  vital  for  ship  operation  are  made  more  reliable  by 
powering  them  from  two  different  sources  through  merging  branches  of  the  pyramid. 
Power  panels  typically  have  multiple  switchboard  inputs.  As  a  result,  isolating  power  to 
a  local  area  in  a  ship  is  unnecessarily  complex.  Standardization  of  the  electric  power 
distribution  system  among  all  ship  classes  is  difficult  at  best  due  to  the  different  number 
of  power  sources  and  vital  loads. 

B.  ZONAL  DISTRIBUTION 

The  zonal  electric  power  distribution  concept  implements  two  power  buses  that 
run  lengthwise  through  the  ship  in  a  configuration  that,  like  the  shipboard  firemain 
system,  maximizes  reliability  and  minimizes  the  possibility  that  both  buses  are  affected 
by  a  single  battle  casualty.  All  operating  generators  are  connected  to  both  buses  via 
switchboards.  Individual  zones  are  designated  typically  to  be  areas  which  directly  relate 
to  flooding  control  zones.  Localized  power  isolation  would,  fiierefore,  be  very  clear-cut 
and  easily  learned  during  a  sailor’s  shipboard  indoctrination  period.  Operation  of  a  zonal 
distribution  system  would  also  be  similar  between  ship  classes. 

C.  DC  ZONAL  DISTRIBUTION 

DC  ZEDS  expands  upon  the  zonal  distribution  concept  in  the  following  manner. 
High  voltage  (approximately  1 100  V)  dc  is  distributed  along  the  main  power  busses.  The 
Ship’s  Service  Converter  Modules  (SSCM)  are  dc-dc  converters  which  buffer  the  main 
bus  and  interzone  loads  as  well  as  regulate  output  voltages.  Where  ac  power  is  required, 
Ship’s  Service  Inverter  Modules  (SSIM)  are  used  to  convert  the  dc  to  ac  and  regulate  die 
ac  voltage  and  frequency  within  the  zone.  Vital  loads  are  connected  to  both  the  port  and 
starboard  SSCM  via  a  diode  auctioneering  process.  Additional  dc-dc  converters 
synthesize  the  regulated  lower  dc  voltages,  for  example  155  V  and  270  V,  required  by 
combat  systems. 
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The  advantages  of  a  DC  ZED  system  are  numerous.  The  most  significant 
advantage  is  increased  survivability  during  casualty  situations.  SSCMs  can  be 
programmed  to  act  as  “smart”  circuit  breakers.  This  would  result  in  maximizing  the 
ship’s  capabilities  to  continue  operating  available  equipment.  A  study  has  shown  that  a 
second  improvement  would  be  weight  and  space  savings  from  the  use  of  solid-state 
components  and  a  reduced  amount  of  cabling  throughout  the  ship  [1].  The  use  of  a  dc 
system  simplifies  operations  such  as  paralleling  loads.  Only  voltages  need  to  be  matched 
on  a  dc  system,  whereas  voltages,  frequency,  and  phase  need  to  be  matched  on  an  ac 
system.  A  fourth  benefit  is  that  variable  speed  ac  motor  control  is  easily  realizable  with 
the  Ship’s  Service  Inverter  Modules  (SSIM),  which  are  programmable  solid-state  dc-ac 
inverters.  Finally,  because  the  operating  generators’  frequencies  are  decoupled  from  the 
ac  loads’  operating  frequencies,  the  prime  movers  can  function  at  their  most  efficient 
speeds. 

Some  disadvantages  to  DC  ZEDS  must  be  acknowledged.  The  state-of-the-art 
solid-state  switches  used  in  power  electronics  today  still  limit  the  dc  bus  voltage  to 
between  1500  volts  and  2000  volts.  Maximizing  bus  voltage  is  important  in  the  effort  to 
minimize  conduction  losses  and  therefore  minimize  conductor  size.  However,  higher 
voltages  pose  concerns  regarding  system  grounding  and  proper  isolation  of  the  dc-dc 
converters.  An  additional  disadvantage  is  that  designers  need  to  more  intelligently  merge 
the  various  components  within  the  system  due  to  stability  concerns. 

D.  PREVIOUS  DC  ZEDS  RESEARCH 

Several  theses  from  the  Naval  Postgraduate  School  have  investigated  various 
aspects  of  DC  ZEDS.  These  areas  include  modeling  portions  of  the  system  to  observe 
controllability  and  stability  issues  of  both  dc-dc  converters  and  synchronous  generators, 
as  well  as  building  and  testing  dc-dc  converters  and  controllers  [2-14].  Most  recently,  an 
investigation  of  frequency-based  load  sharing  between  current-mode-controlled  buck 
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converters  was  conducted  which  included  an  rms  frequency  estimation  algorithm  and 
accompanying  hardware  to  achieve  the  desired  load  sharing  characteristics  [15]. 

Previous  work  conducted  in  the  fields  of  synchronous  machine  modeling  and  dc 
link  filter  stability  issues  were  used  as  the  foundation  for  this  thesis.  Modeling  of  an 
uncontrolled  synchronous  machine  in  a  voltage-behind-reactance  representation  was 
conducted  by  Sudhoff  and  Wasynczuk  [16]  and  Branson,  Ciezki,  Sudhoff,  and 
Wasynczuk  [17].  Research  conducted  by  Sudhoff,  Corzine,  and  Glover  at  die  University 
of  Missouri-Rolla,  and  Hegner  and  Robey  at  the  Naval  Surface  Warfare  Center  makes 
observations  on  stability  conditions  in  a  simplified  dc  link  model  [18].  A  technical  report 
by  Ciezki  and  Ashton  at  the  Naval  Postgraduate  School  discussed  converter  and  converter 
control  design  methodology  [19].  A  dissertation  by  Belkhayat  explored  an  impedance- 
based  stability  criteria  for  a  stand-alone  dc-dc  converter  [20].  This  small-signal 
requirement  was  investigated  by  application  of  Nyquist  stability  criteria.  Finally, 
research  conducted  by  Belkhayat  and  Cooley  at  the  Naval  Surface  Warfare  Center,  and 
Witulski  at  the  University  of  Arizona  investigated  large-signal  stability  criteria  for 
distributed  systems  -with  constant  power  loads  [21]. 

E.  THESIS  GOALS 

This  thesis  explores  various  stability  criteria  for  a  representative  ship-board  dc 
power  distribution  system  that  includes  a  synchronous  generator,  a  phase-controlled 
rectifier,  a  LC  Imk  filter,  and  a  dc-dc  converter.  Specific  limits  on  the  filter  parameters 
and  converter  input  impedance  which  can  be  given  as  design  specifications  to  contractors 
are  the  desired  results.  The  analysis  methods  examined  are  a  simplified  constant-power 
dc  link  model,  a  model  based  on  system  impedances,  and  a  Lyapunov-based,  large-signal 
nonlinear  representation.  All  models  were  studied  using  the  Advanced  Continuous 
Simulation  Language  (ACSL)  together  with  specially  written  MATLAB  script  files.  In 
addition,  an  ACSL  model  of  the  entire  system  was  created  for  more  advanced  analysis. 
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F. 


CHAPTER  OVERVIEW 


The  derivation  of  modeling  equations  for  a  reduced-order  synchronous 
machine/rectifier  drive  are  detailed  in  Chapter  11.  The  accompanying  model  code  is  then 
explained  in  detail.  Chapter  III  examines  a  simplified  dc  link  model  and  the  limitations 
placed  to  ensure  stability  by  the  Hurwitz  criterion.  The  results  of  simulations  used  to 
verify  the  limitations  are  presented.  Chapter  IV  then  examines  stability  criteria  based  on 
a  system  impedance  model.  Design  of  a  dc  converter  is  outlined,  the  required  impedance 
values  are  determined  and  the  Nyquist  criterion  is  applied  to  the  model.  Nyquist  plots  for 
stable  and  unstable  conditions  are  displayed.  Chapter  V  covers  the  application  of  a 
Lyapunov-based  large-signal  stability  method  to  a  simplified  system  model  to  ensure  that 
large  system  disturbances  are  eveiluated  in  addition  to  the  earlier  small-signal  models. 
Specifically,  the  equilibrium  point  criterion  and  mixed  potential  criterion  are  evaluated. 
Finally,  Chapter  VI  presents  the  conclusions  and  discusses  areas  available  for  further 
research. 
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II.  SYNCHRONOUS  MACHINE  -  RECTIFIER  MODELING 


A.  REDUCED-ORDER  MACHINE  DEVELOPMENT 

The  first  step  in  developing  an  overall  model  is  accurately  modeling  the  electrical 
source,  which  in  the  system  being  evaluated  is  a  three-phase,  ac  synchronous  generator 
connected  to  a  phase-controlled  rectifier.  A  detailed,  valve-by-valve  simulation  is 
computationally  intensive  and  results  in  considerably  long  run  times.  Since  the  dynamic 
effects  of  the  system  are  important  for  the  analysis,  a  time-efficient  reduced-order  model 
is  used  for  simulation  and  analytical  purposes.  Figure  II-l  shows  a  synchronous  machine 
depicted  as  a  voltage-in  /  current-out  model.  The  machine  terminal  voltages,  therefore, 
are  required  as  both  inputs  to  the  generator  governing  equations  and  for  use  in 
determining  the  commutation  properties  of  the  thyristors  in  the  rectifier.  Thus,  there  is  an 
inherent  algebraic  incompatibility  that  must  be  resolved  by  the  model  synthesis. 

1.  Approach 

The  following  methodical  approach  was  used  in  the  development  of  the  reduced- 
order  model.  The  machine  stator  electrical  transients  (known  as  the  “fast”  transients) 
were  ignored  for  simulation  purposes.  A  voltage-behind-reactance  (VBR)  representation 
of  the  synchronous  machine  was  then  derived.  The  VBR  model  was  subsequently  used 
.along  with  the  different  system  modes  to  evolve  an  equation  for  the  average  direct  current 
output  voltage  of  the  rectifier  (V<,c).  The  source  representation  was  completed  by 
establishing  the  required  synchronous  machine  stator  currents. 
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2.  Synchronous  Machine  Modeling  Equations 


From  reference  [16],  the  equations  used  to  describe  a  reduced-order  synchronous 
machine  are  given.  These  equations  were  derived  by  setting  the  derivatives  of  the  stator 
flux  Imkages  per  second  equal  to  zero  and  fixing  the  reference  frame  speed  equal  to 
(the  electrical  angular  velocity).  Further,  assuming  an  isolated  synchronous  machine  in 
this  case,  we  find  that  0)^  —  cOj.  (the  rotor  electrical  angular  velocity).  The  stator  voltage 
equations  are  as  follows: 


v''  = 

''qs 


r  c  +  —  'f: 


s  "qs 


v;  =  4  + 


0)^ 


ds 


(1) 

(2) 


where  v',  and  are  the  machine  terminal  voltages  expressed  in  the  rotor  reference 
frame  (RRF),  and  are  the  stator  currents  in  the  RRF,  and  are  the  stator 
fliix  linkages  per  second  in  the  RRF,  and  co^  is  the  base  electrical  angular  velocity  of  the 
machine.  The  transformation  matrix  which  relates  the  RRF  quantities  to  the  abc 
reference  frame  quantities  is: 
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where  r^,  and  r  'jy  are  the  referred  winding  resistances  of  the  rotor  coils.  The 

referred  values,  designated  by  the  primed  quantities,  are  simply  scaled  versions  of  actual 
machine  quantities.  The  scaling  factor  is  related  to  the  ratio  of  winding  turns  between  the 


rotor  and  the  stator  (see  [22].). 

Stator  flux  linkage  per  second  equations  are  given  by: 

%  =  +  (7) 

Rotor-referred  flux  linkage  per  second  equations  are  expressed  by: 

=  +  (9) 

~  ^Ikd  ^kd  ^ md  ^ds  "^^kd  (1^) 

^fd  ~  ^Ifd  4  ^  ^md  ^ds  '^^kd  (1^) 
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where  the  various  reactances  are  constants  for  a  given  machine  provided  operation 
remains  magnetically  linear. 

One  expression  for  the  torque  may  be  written  in  terms  of  the  stator  currents  and 
stator  flux  linkages  per  second: 

(12) 

where  P  is  the  number  of  poles  of  the  machine. 


The  state  variables  for  the  system  are  ,  o, ,  and  .  The  stator  flux 

linkages  per  second  are  not  state  variables  since  we  have  ignored  the  derivative  terms  in 
those  voltage  equations.  The  next  step  is  to  express  and  in  terms  of  state 

variables  and  stator  current.  Solving  (10)  and  (1 1)  for  4  and  4  and  solving  (9)  for  4 . 
we  find: 


1 

1 

- \ 

1 

1 

p- 

(13) 


/• 

'*9  Y' 
^kq 


X, 


X 


mg  -r 
*gs 
kq 


(14) 


Equations  (13)  and  (14)  are  substituted  into  the  stator  flux  linkage  per  second 
Equations  (7)  and  (8)  to  yield: 


(15) 

(16) 

where  the  q  and  d-axis  subtransient  reactances  are  denoted  and  Xj ,  and  the  q  and  d- 
axis  subtransient  flux  linkages  per  second  are  denoted  'Pj  and  'Pj .  It  can  be  noted  that 

the  subtransient  flux  linkages  per  second  are  functions  of  the  state  variables  of  the 
synchronous  machine.  The  next  step  is  to  then  substitute  Equations  (15)  and  (16)  into  the 
stator  voltage  Equations  (1)  and  (2)  to  give: 
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'^qs  ’  s^qs 


^ds  -  ^s^ds 


.±r 

CO. 


G). 


It  is  customary  to  define: 


■-xX+^% 

(17) 

b 

(18) 

(19) 

1 

k 

(20) 

^4  . 

then  Equations  (17)  and  (1 8)  can  be  rewritten: 


v:.=rM-^XJl+E: 


(21) 

(22) 


Equations  (21)  and  (22)  are  the  desired  reduced-order  voltage-behind-reactance  (VBR) 
representation  of  the  synchronous  machine. 


3.  Modeling  Average  DC  Voltage 

The  next  step  involved  in  modeling  the  electrical  power  source  is  to  derive  an 
expression  for  the  average  dc  voltage  using  the  VBR  model  together  with  boundary 
values  associated  with  the  conduction  of  the  rectifier.  In  reference[16],  the  technique 
used  was  to  find  the  average  dc  voltage  during  a  single  switching  interval,  since  each  of 
the  six  steady-state  switching  intervals  in  a  three-phase  rectifier  are  identical.  During  a 
single  switching  interval,  the  average  dc  voltage  is  a  function  of  the  phase  voltages  of  two 
of  the  phases  of  the  synchronous  machine  over  the  switching  interval  as  shown  by: 

(23) 

/  3 
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where  p  is  the  firing  delay  angle  from  the  natural  point  of  commutation  with  respect  to 
the  rotor  position  [16]. 

The  VBR  model  can  be  used  to  show  that: 


3  CO 

7t  o*  er=T*P 


(24) 


The  small  effect  of  the  stator  resistance  is  ignored. 

The  flux  linkages  per  second  can  then  be  determined  at  the  two  limits  using  the 
DC  cxirrents  at  each  limit  to  yield  the  following  equation  for  the  average  dc  voltage: 


CO. 


3  CO. 


Tt  (O^ 


cosa-- 


V  o 


(25) 


where  a  is  the  firing  delay  angle  from  the  natural  point  of  commutation  with  respect  to 
the  source  electrical  angle  [16]. 


4.  Calculation  of  Stator  Currents 

The  final  step  in  modeling  the  power  source  is  to  calculate  the  stator  currents. 
This  is  necessary  since  the  referred  rotor  currents  (Equations  (13)  and  (14))  are  expressed 
in  terms  of  rotor-referred  flux  linkages  per  second  and  rotor  reference  frame  stator 
currents.  In  the  developed  model,  the  rotor  reference  frame  stator  currents  are  needed  to 
update  the  state  of  the  machine  during  each  time  step  iteration.  The  non-ideal  case  of  the 
current  flow  through  the  six  individual  thyristors  as  a  function  of  the  three-phase 
synchronous  generator  output  as  shown  in  Figure  II-2.  As  it  would  be  expected,  these 
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1  ^ 


/3  f*^ 

The  non-ideal  switching  interval  can  further  be  divided  into  two  distinct  events.  The  first 
part  of  the  switching  interval,  called  the  conduction  period,  consists  of  the  section  of  time 
that  there  is  only  one  device  conducting  current  on  each  rail.  The  second  part  of  the 
switching  interval,  called  the  commutation  period,  occurs  when  current  is  being 
transferred  from  one  device  to  another  device  on  the  same  rail.  Thus: 

^qs,ave  ^qs,com  ^qs,cond  (28) 

^ds,ave  ~  ^ds,com  ^ds,cond  (29) 

If  u  is  defined  as  the  conduction  angle,  then  Equations  (26)  and  (27)  can  be 
subdivided  into  the  conduction  and  commutation  periods  as  follows: 


iqs,com=-^  \ils(6r)d0r 


3 


"i  x-f 


1  ~ 

iqs,co„,=^  \iqs(0r)de. 


^3  f 


j  T  ^ 

^*,conrf  =  P*  (.0r)d0r 


3  f^P^u 


Since  the  currents  are  known  at  the  boundaries  of  the  conduction  interval, 
modeling  equations  for  the  conduction  currents  are  easily  obtained  from  Equations  (32) 
and  (33): 
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The  only  unknown  in  the  above  equations  is  the  commutation  angle  (u).  This  can 
be  found  while  examining  the  commutation  period  dynamics.  Finding  the  commutation 
currents  is  a  more  difficult  endeavor.  During  commutation,  two  of  the  phases  are  shorted 
together.  As  a  result,  the  difference  between  the  flux  linkages  per  second  of  the  two 
phases  is  a  constant,  specified  in  [16]  as  A", .  The  quantity  can  be  solved  for  in  the 
abc  reference  firame  yielding: 


cosf  e,  +  %  + 1 J 

(x;  -  x;  )cosf  + (x;  +  x; ) 

(x;  -  x;  )cos[^2e, + + i  (x;  +  x; ) 

,(x;  -x;)co/2e, +1  ]+fe  +x;) 

where  the  aforementioned  is: 

K,  =s(-Xmsfi+'¥:smfi]  (x;  -x;)cosf2^+yj-i(x; +x;)  /* 

(37) 

When  commutation  is  complete,  as  a  function  of  0^  is  zero.  An  implicit 
expression  for  the  commutation  angle  u  is  found  by  setting  equal  to  zero  in  Equation 

(36): 

-V3[-'f;smC8  +  !/)+T;cosC8  +  «)] 

-  (X,  -Xjcos[^2;5  +  2«-^J  +  i(x,  +Xj 
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Given  the  state  of  the  machine,  the  dc  link  and  the  firing  angle,  the  commutation 
angle  can  be  detennined  by  setting  Equation  (38)  equal  to  zero  and  solving  for  u .  This  is 
done  in  the  simulation  by  iteration  to  achieve  a  commutation  angle  within  0.001  radians. 
The  iterative  method  used  simply  involves  halving  the  range  in  which  the  commutation 
angle  exists  during  each  step  until  the  range  is  narrowed  to  the  above  stated  tolerance. 
Once  u  is  calculated,  the  conduction  current  components  can  be  determined  using 
Equations  (34)  and  (35). 

The  q  and  d-stator  currents  in  the  RRF  are  defined  during  the  commutation  period 


as; 


4  = 


2V3 

3 

2V3 


r 


»a.COS 


*a.sm 


9.+- 


+  4,cos6>. 


(39) 


(40) 


V 

Eqioatipn  (36)  is  substituted  into  the  above  equations,  which  are  then  substituted 
into  Equations  (30)  and  (31).  The  resulting  integrals  cannot  be  solved  by  analytical 
methods.  Instead,  a  numerical  method  of  integration  called  the  Simpson  method  is  used 
[16].  The  Simpson  method  yields  the  following  equations  for  the  commutation  currents 
over  the  commutation  angle  by  utilizing  a  weighted  measurement  of  the  q  and  d-currents 
at  increments  through  the  commutation  process,  including  the  boimdaries: 


qSyCom 


^ds^com 


3  u 

^12 


3  u 
7t  12 


'qs 


+4/1 


(41) 


'■  3 


+- 
3  4 


3  2  3  2 


(42) 

Now  that  both  the  conduction  and  commutation  currents  are  known,  the  average 
RRF  stator  currents  are  calculated  with  Equations  (28)  and  (29).  These  currents  are  then 
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injected  into  the  synchronous  machine  rotor  flux  linkage  per  second  equations  shown  in 
Equations  (15)  and  (16). 

B.  DESCRIPTION  OF  ACSL  CODE  AND  SIMULATION 

The  derived  model  is  programmed  in  the  Advanced  Continuous  Simulation 
Language  (ACSL)  .  A  listing  of  the  .csl  and  .cmd  files  is  documented  in  Appendix  A. 
The  .csl  file  contains  the  equation  description  of  the  synchronous  generator-phase 
controlled  rectifier  in  ACSL  statement  format.  The  .cmd  file  contains  sets  of  runtime 
instructions  to  manipulate  the  model  output  into  a  usable  form.  ACSL  was  chosen  for  its 
execution  speed  and  its  strength  and  ease  of  working  with  and  sorting  continuous 
modeling  equations.  The  feature  is  desirable  due  to  the  number  of  variables  involved 
with  the  model  as  well  as  the  interconnectivity  of  the  numerous  variables. 

1.  .CSL  File  Description 

A  block  diagram  representation  of  the  model  in  terms  of  the  required  inputs  and 
resultant  ouqjuts  is  shown  in  Figure  II-3.  The  initio  section  of  the  .csl  file  holds  the 
machine  constants  used  in  the  test  model.  The  assumed  per  unit  electrical  machine 
parameters  were  taken  from  reference  [17]  for  the  following:  ,  X,^ ,  ,  X^ ,  , 

X]^ ,  X'lj-j ,  X',^ ,  r'^ ,  r'jj ,  r'^  and  ,  where  the  leakage  reactance  terms  with  a 
subscript  are  used  to  calculate  the  q  and  d-axis  subtransient  reactances  {X^  and  Xj). 

The  derivative  section  contains  the  equations  that  describe  the  system.  Within  the 
derivative  section,  there  are  two  subroutine  calls  and  two  function  calls.  The  two 
FORTRAN  subroutines  implement  the  Simpson  method  described  in  reference  [16]  to 
calculate  the  five  different  current  components  in  equations  (40)  and  (41).  The  first 
FORTRAN  function  calculates  the  commutation  angle  (m  )  as  described  in  an  earlier 
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Figure  11-3  Block  Diagram  of  Synchronous  Generator-Phase  Controlled  Rectifier 

ACSL  Model 

section.  The  second  FORTRAN  function  calculates  the  firing  delay  angle  widi  respect  to 
rotor  angle  {/3)  using  the  same  method  as  adopted  for  calculating  u . 

The  model  accoimts  for  the  fact  that  p  must  be  set  such  that  the  thyristor  is  not 
gated  until  after  it  is  forward  biased  by  constraining  p  to  a  minimum  value,  ,  that 
denotes  the  earliest  that  a  thyristor  can  be  turned  on.  The  firing  angle  for  the  thyristor 
with  respect  to  the  point  its  anode-to-cathode  voltage  becomes  positive  is  then 
represented  by  a ,  which  is  the  difference  between  P  and  . 

2.  .CMD  File  Description 

The  .cmd  file  is  configured  to  run  a  simulation  to  produce  data  points  which 
define  and  as  a  function  of  p .  The  resulting  graphs  for  1^^  =  0.9  pu  are  shown  in 
Figure  II-4.  The  effect  of  this  is  constant  values  for  and  at  values  of  p  less  than 
the  minimum  value  that  represents  gating  the  thyristor  prior  to  it  being  forward  biased.  It 
can  be  seen  firom  the  graphs  that  p^-^  is  at  approximately  30°.  With  p  equal  to  p^:^ , 
and  are  at  their  extreme  per  unit  values  of  1 .4  and  -0.85  respectively.  and 
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change  linearly  as  is  increased  in  value  when  it  is  greater  than  decreases 

from  a  positive  value  and  increasing  from  a  negative  value,  corresponding  to  generator 
action.  When  exceeds  approximately  70°,  the  average  dc  voltage  becomes  negative 
and  the  average  torque  becomes  positive,  corresponding  to  motor  action.  This  linear 
region  lasts  until  P  reaches  approximately  105°. 


Figure  II-4  Beta  vs  Vdc  and  Te 

A  detailed  description  of  a  model  for  the  system  source  has  now  been  presented. 
The  next  chapter  discusses  the  evaluation  of  a  simplified  dc  link-converter  model  for 
stability  criteria  in  terms  of  major  system  parameters. 
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III.  SIMPLIFIED  DC  LINK/CONVERTER  MODEL 


A.  DEVELOPMENT  OF  THE  MODEL 

One  characteristic  predominant  in  dc-dc  converters  is  the  apparent  small-signal 
negative  input  impedance  which  is  a  result  of  the  high  bandwidth  regulation  that  makes 
the  converter  appear  as  a  constant  power  sink.  In  Ref.  [18],  it  is  noted  that  one  method 
that  can  be  used  to  ensure  stability  is  to  increase  the  dc  link  capacitance.  This  is  not  a 
desirable  remedy  since  it  requires  the  use  of  more  space  and  increases  the  weight 
contribution  of  the  overall  distribution  system.  In  addition,  since  electrolytic  capacitors 
have  relatively  low  reliability,  large  banks  are  unattractive  owing  to  the  time-intensive 
task  of  identifying  a  shorted  capacitor.  The  authors  in  [18]  develop  a  simplified  model, 
as  shown  in  Figure  III-l,  to  illustrate  the  cause  of  negative  input  impedance  instability. 


Figure  m-l  Simplified  DC  Link  Model 


Several  simplifying  assumptions  are  made.  First,  the  prime  mover  and  voltage 
regulator  dynamics  are  assumed  to  have  significantly  long  time  constants.  Second,  all 
losses  in  the  machine  and  rectifier  are  assumed  to  be  zero.  Third,  the  rectifier  is  assumed 
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to  be  uncontrolled  and  finally,  field  voltage  regulation  maintains  the  source  voltage  equal 
to  the  commanded  source  voltage  resulting  in  no  need  for  subtransient  information. 

The  equations  which  describe  the  non-linear  simplified  model  are  then  found  to 
be: 


dt 


L 


=  -—^da—rh  + 


L 


L 


dv 


p* 


dci  _  ^  j _ 

dt  Cg  ^e^dci 


(43) 


(44) 


where  is  a  function  of  the  filter  inductance  and  the  transient  commutating  inductance; 

is  a  function  of  the  commutating  inductance,  machine  stator  resistance,  and  link  filter 

resistance;  and  includes  the  filter  capacitance  and  any  input  capacitance  associated 

with  the  converter.  P  is  the  assumed  constant  power  being  drawn  by  the  converter. 
Linearization  of  the  above  equations  result  in: 


d  Ai, 


1 


R. 


^  =  -—Av,^,—^Ai,+  —  Av. 


dt  L, 
d  Av . .. 


1  A-  ^  A 

2  ^^dci 

Ce  C,V^i 


(45) 

(46) 


Furthermore,  the  linearized  converter  input  current  is: 


^4/  = 


P* 


dci  L  *^dci  J 


-P* 


Av 


dci 


(47) 


''dci 


Thus,  the  small-signal  converter  input  impedance  may  be  found  from: 


Z. 


^dci 

^idci~  P* 


(48) 


Therefore,  for  a  given  power  being  supplied  by  the  converter  and  a  given  input  voltage 
,  the  input  impedance  must  be  negative.  It  is  this  negative  resistance  which  sets  up 
the  possibility  for  oscillations  with  the  input  LC  filter. 
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B.  EVALUATION  OF  THE  MODEL 

Confinnation  that  the  negative  input  impedance  results  in  system  instability  is 
obtained  by  examining  the  characteristic  equation  of  the  model.  In  order  for  the  system 
to  be  stable,  a  necessary  and  sufficient  condition  of  the  Hurwitz  criterion  requires  that  all 
coefficients  of  a  second-order  characteristic  equation  be  positive.  From  Equations  (45) 
and  (46),  the  characteristic  equation  for  the  system  is: 


4  Qvi, 


4Q  ^e^e^dci 


The  resulting  stability  requirements  are  then: 


4  P* 

—  > - r 

L.  Cvi 


^>4  (51) 

Equation  (51)  is  typically  satisfied  by  matching  the  load  characteristics  to  the  given 
source  resistance.  Thus,  Equation  (50)  becomes  a  first-order  design  criterion  for  the  link 
filter  parameters.  A  rearranged  version  of  Eqiiation  (50)  illustrates  that  we  can  ensure 
small-signal  stability  by  including  additional  bus  capacitance  such  that: 

C  P* 

—  >-V  (52) 

4  4vL 

ACSL  code  for  the  .CSL  and  .CMD  files  which  model  both  the  linearized  and 
non-linearized  representations  of  the  model  is  given  in  Appendix  C.  Simulations  were 
conducted  in  which  P*  was  assumed  to  be  200  kW,  was  assumed  to  be  1  lOOV, 
was  assumed  to  be  0.01  Q ,  and  was  assumed  to  be  0.5  mH.  This  resulted  in  requiring 
Q  to  be  8265  //  F  to  just  be  stable  in  the  event  of  large  disturbances. 

The  simulated  output  of  a  transient  from  10%  rated  power  to  100%  rated  power  at 
time  0.5  seconds  and  a  return  to  10%  rated  power  at  time  3  seconds  is  shown  in  Figure 
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III-2.  A  value  of  C^=  8500  juFwas  used  to  better  show  the  stable  response  due  to  the 
low  damping  of  this  configuration.  The  system  can  be  seen  to  be  stable  after  both  the  up- 
power  and  down-power  transients.  The  simulated  output  of  a  similar  transient  sequence 
with  the  down-power  transient  occuring  at  time  0.6  seconds  is  shown  in  Figure  III-3.  For 
this  simulation,  was  lowered  to  826.5  /liF  .  It  is  shown  that  following  the  up-power 

transient,  the  system  becomes  unstable.  Stability  is  not  recovered  following  the  down- 
power  transient. 

We  have  now  established  the  governing  stability  criteria  dictated  by  a  simplified 
dc  Imk/converter  model.  In  the  next  chapter,  we  examine  a  more  detailed  model  that 
defines  the  requirements  for  stability  using  subsystem  impedances. 


Figure  111-2  Up-Power,  Down-Power  Transient,  C.  =  8500  MicroFarads 
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IV.  IMPEDANCE-BASED  STABILITY  CONSTRAINTS 


A.  STABILITY  CRITERIA  USING  SUBSYSTEM  IMPEDANCES 

One  analysis  that  can  be  conducted  to  determine  the  stability  of  the  overall  system 
involves  determining  the  input  impedance  of  the  dc  converter  and  the  output  impedance 
of  the  filter.  Together,  the  poles  of  the  overall  system  transfer  fimction  are  foimd.  The 
transfer  fimction  can  then  be  examined  to  determine  if  the  Nyquist  stability  criteria  are 
met.  The  steps  involved  in  this  process  are:  design  the  converter  to  be  used,  determine 
the  converter  controller  gains,  calculate  the  converter  input  impedance,  calculate  the  filter 
output  impedance,  and  use  the  results  to  analyze  stability.  For  analysis  purpose,  two 
converters  are  developed  for  use  in  simulations.  The  first  is  designed  for  a  low  switching 
firequency  used  for  a  high-power  hard  switching  process.  The  second  will  be  designed  for 
a  higher  switching  fi-equency  used  for  a  high-power  soft  switching  process. 

B.  DESIGN  OF  THE  DC  CONVERTER 

The  layout  of  a  common  step-down  dc-dc  converter  in  a  Buck  Chopper 
configuration  is  shown  in  Figure  IV-l.  To  design  the  converter,  the  first  issue  which 
must  be  resolved  is  the  system  requirements.  To  establish  a  baseline  for  the  design,  we 
require  maintaining  continuous  current  conduction  mode  for  all  loads  that  are  greater  than 
ten  percent  of  the  rated  power.  For  this  thesis,  the  design  specifications  listed  in  Table 
rV-1  are  assumed. 
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Specification 

Value 

Rated  Power 

200  kW 

Filter  Output  Voltage  (V^j.) 

1100  V 

Load  Voltage  (V„J 

950  V 

Ripple  Voltage  (A  V„J 

5  V 

Hard  Switching  Frequency  (fshard) 

4  kHz 

Soft  Switching  Frequency  (fj  soft) 

20  kHz 

Table  IV-1  Design  Specifications  for  the  DC-DC  Converter 


R  (load) 


Figure  IV-1  DC-DC  Converter  -  Buck  Chopper 
1.  Design  of  the  Components 


The  authors  in  [19]  outline  a  sound  methodology  for  converter  design.  The 
nominal  duty  cycle  (DJ  of  the  switch  is  defined  as  the  percentage  of  time  that  the  switch 
is  closed.  For  the  designed  converter,  assuming  continuous  current  conduction,  we  get: 


I>„  = 


950F 

HOOF 


0.8636 


(53) 
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Next,  based  on  the  rated  power  and  the  desired  output  voltage,  the  rated  load  resistance 
can  be  found  by: 


R 


rat 


VL  (950ry 

P^,  200/cW 


4.5125Q 


(54) 


The  critical  resistance  (R^„; )  is  the  load  resistance  at  the  boundary  between  continuous 
and  discontinuous  current  conduction  modes.  For  this  converter: 

R^„.,  =  lOx  =  (10)x(4.5125Q)=  45.125Q  (55) 

The  main  components  of  the  converter  can  now  be  selected.  The  hard  switching 
converter  parameters  are  determined  first.  The  critical  inductance  (1^^,)  is  the  minimum 
inductance  required  to  maintain  continuous  current  conduction  mode  during  all 
conditions  required  by  the  design  specifications.  The  critical  inductance  is  found  by: 

(56) 

In  order  to  maintain  a  comfortable  margin,  L  is  assumed  to  be  ImH. 

The  minimum  capacitance  is  designed  around  ensuring  that  the  steady-state 
output  voltage  ripple  specification  is  met: 

fc-v»)  =  202.4/J’  (57) 

® ^crit  J  Syhard out 

Once  again,  a  margin  is  included  to  ensure  that  the  capacitor  is  large  enough  to  provide 
sufficient  energy  during  a  transient  to  maintain  output  voltage.  The  capacitor  selected 
should  not  be  too  large,  however,  because  this  would  result  in  small  changes  in  output 
voltage  causing  the  duty  cycle  to  oscillate  between  full-on  and  full-off.  C  is  assumed  to 
be  2mF. 


2.  Converter  Controller  Design 

The  buck  chopper  exhibits  an  oscillatory  response  as  a  result  of  the  undamped 
L-C  filter  at  the  converter  output.  To  ensure  that  the  converter  operates  satisfactorily 
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during  transients  and  to  provide  steady-state  converter  regulation,  a  multi-loop 
feedforward  controller  is  utilized  for  duty  cycle  control.  The  governing  controller 
equation  is  [19]: 

d{>)  =  -^  +  KK,{r^  -  v„ )+  h.K,  j(v„j  -  v„, )dt  +  ) 


(58) 

where  h^,  h„,  and  h^  are  the  feedback  gains,  and  K;  are  the  voltage  and  current  sensor 
gains,  V„f  is  the  desired  output  voltage,  and  is  the  sensed  output  current.  This 
controller  design  assumes  that  the  converter  is  operating  autonomously.  Additional  terms 
are  required  to  parallel  multiple  dc-dc  converters.  The  first  component  of  Equation  (58) 
provides  feedforward  action  to  ensure  rapid  response  to  disturbances  in  the  input  voltage. 
The  second  and  third  terms  implement  a  proportional-plus-integral  (PI)  outer  voltage 
control  while  the  fourth  term  realizes  an  inner  current  control.  Effectively,  this  last 
component  represents  the  feedback  of  the  derivative  of  the  capacitor  voltage  so  that  a 
PID-t3rpe  contol  is  realized. 

By  assuming  that  the  desired  output  voltage  )  remains  constant  and  the 
converter  load  is  purely  resistive  in  nature,  we  find  that  the  operation  of  the  converter  can 
be  described  by  a  set  of  averaged  linear  ordinary  differential  equations  with  the  following 
characteristic  polynomial: 


+d^s^  +d,s  +  do  (59) 

where  dj,  dj,  and  dg  (fimctions  of  the  converter  parameters  and  control  gains)  describe  the 
system  pole  locations  and,  ultimately,  the  converter’s  transient  response.  Upon  selecting 
desired  characteristic  polynomial  coefficients,  we  establish  gain  relationships  as  follows: 


^dc 


d,- 


LC. 


KK,.  = 


d^LC 


dc 


(60) 

(61) 
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(62) 


h,K,  = 


-L 


V, 


dc 


d.- 


RC, 


The  utilization  of  a  third-order  Bessel  function  polynomial  as  a  system  prototype 
is  known  to  result  in  a  desirable  system  response.  In  addition,  it  minimizes  the  distance 
that  the  open-loop  poles  must  be  moved  and  hence  minimizes  the  control  effort.  The 
specific  poles  for  the  prototype  system  are  located  at: 

^1,2  =  -0.7455®  y0.71 12®  (63) 


=  -0.9420® 


(64) 


where  ®  is  the  desired  closed-loop  bandwidth.  The  radian  frequency  (® )  must  meet 
three  criteria  to  provide  optimized  results.  The  first  criteria  is  that  ®  must  be  at  least  one 
decade  below  the  switching  frequency  to  prevent  undesirable  controller-switch 
interactions.  The  second  criteria  is  that  ®  should  be  minimized  to  prevent  additional 
noise  fi:om  being  introduced  into  the  system  through  the  controller’s  operation.  The  final 
criteria  is  that  ®  must  be  large  enough  to  achieve  the  desired  settling  time  following  a 
transient.  As  a  result  of  these  criteria,  ®  was  chosen  for  the  hard  switching  converter  to 
be  1500  ”^/sec.  The  poles  for  the  designed  converter  are  then: 

.y,,  =-1118.25^71066.8 


^3  =-1413 


The  system  characteristic  equation  is  then: 

^'+3639.55^  +5526354.85  +  3351129059.4 
and  the  gain  selections  become: 


h,K,  =  0.00914 


h„K,  =6.093 


(65) 

(66) 

(67) 

(68) 

(69) 

(70) 


=  -0.003299 

Similar  calculations  were  made  for  a  soft  switching  buck  chopper  with  a 
switching  frequency  of  20  kHz.  was  1 53.9  //  H,  with  a  value  of  200  //  H  chosen  to 
provide  a  margin.  was  found  to  be  40.5  //  F,  with  a  value  of  400  ju  F  chosen  to 
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ensure  adequate  energy  storage  capacity.  The  desired  system  poles  utilizing  the  Bessel 
prototype  with  co  =  3000  ”‘*4^  are: 


^,2  =-2236.5  2133. 6  (71) 

^3  =  -2826  (72) 

The  resultant  gains  are: 

h,K,  =  0.00070508  (73) 

4^,  =1.9636  (74) 

4^:,.  =-0.001317  (75) 


3.  FOter  -  Converter  Model  in  ACSL 

The  .csl  and  .cmd  files  for  an  ACSL  model  of  a  dc  link  filter  attached  to  the  input 
of  a  dc-dc  converter  with  feedback  control  are  contained  in  Appendix  D.  The  model 
specifically  contains  the  component  parameters  derived  for  the  hard  switching  converter. 
It  assumes  an  ideal  dc  input  voltage.  The  initial  conditions  for  the  simulation  are  that  the 
filter  capacitor  is  fully  charged,  and  the  converter  is  secured.  Converter  operation  is 
initiated  by  ramping  up  the  reference  voltage.  Converter  control  is  implemented  with  the 
method  outlined  in  the  previous  section. 

Several  simulations  were  run  with  various  filter  parameters.  The  behavior  of  the 
system  with  changes  in  the  filter  capacitance  were  of  particular  interest.  When  filter 
parameters  were  within  those  specifed  by  the  stability  requirements,  the  output  was  stable 
as  expected.  However,  of  interest  was  the  observation  of  a  region  of  stability  outside  of 
the  criteria  bounds  established  in  the  stability  analysis.  The  band  was  observed  to  exist 
when  filter  capacitor  values  between  approximately  25  and  1 10  jjF  were  used  in  the 
simulations.  It  is  believed  that  this  region  of  stability  is  a  result  of  the  feedback  control 
for  this  particular  set  of  parameters.  In  general,  changes  in  the  filter  capacitance  resulted 
in  changes  in  the  ripple  voltage  magnitude  and  frequency. 
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c. 


IMPEDANCE  CALCULATIONS 


A  representation  of  a  source-filter-converter-load  is  shown  in  Figure  IV-2,  as  also 
shown  in  reference  [20].  The  overall  transfer  function  of  this  system  is  shown  to  be: 


V,  \  +  ZJZ, 

where  is  the  output  impedance  of  a  line  LC  filter  and  source,  and  Z,  is  the  input 
impedance  of  a  loaded  converter  with  closed-loop  control.  Stability  of  the  system  can  be 
determined  by  analyzing  the  poles  of  this  transfer  fimction.  This  section  will  describe 
how  Z^  and  Z,  are  determined  for  this  system  as  well  as  analyze  the  behavior  of  the 
poles  of  the  system. 


Zs  z, 

Figure  IV-2  DC-DC  Converter  With  Input  Filter 


1.  Determination  of  Converter  Input  Impedance 


Due  to  the  non-linear  nature  of  a  buck  chopper,  determination  of  the  input 
impedance  needs  to  be  performed  on  a  linearized  version  of  the  model.  Referring  back  to 
Figure  IV-1,  we  find  the  input  current  is  related  to  the  averaged  inductor  current  by: 

h  =^h  ,  ('77) 

where  d  is  the  converter  duty  cycle.  Linearization  of  Equation  67  yields: 

Ai^=d„Ait^+ij^„Ad  (78) 
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where  bd  is  the  perturbation  variable  given  by  : 

M^d-d„  (79) 

The  constants  d^  and  are  operating  point  vdues  for  those  variables.  By  assuming 
is  zero,  the  linearization  of  Equation  48  can  be  shown  to  yield: 

bd  =  K^bv^  +  r,  (5)Av^  (80) 

where  AT,  is: 


-K 


(81) 


//2,0 


and  T^(js)  is: 


T  ( j')  =  ^  +h^s  +  h„ 


(82) 


Further  definitions  that  are  made  involve  the  transfer  fimction  between  the  variation  in 
input  voltage  ( Av;„ )  and  the  variation  in  the  state  variable  inductor  current  ( A/^ )  as  well 
as  the  variation  in  the  state  variable  capacitor  voltage  ( Av  These  relationships  are: 

A4=7’2(5)Av,„  (83) 

Av,=r3(^)Av,„  (84) 

The  MATLAB  code  used  to  determine  the  coefficients  for  J’2  (j)  and  (5)  is 
contained  in  Appendix  E.  The  method  used  to  find  these  involves  deriving  the  state- 
space  form  of  the  linearized  system  equations,  then  utilizing  the  state-space  to  transfer 
function  (ss2tf)  MATLAB  function  with  a  C  matrix  that  zeros  all  outputs  except  bi^^  for 
the  T^is)  calculation  and  Av^  for  the  r3(5)  calculation. 

The  desired  input  impedance  is  defined  by: 

Av,„ 


Z,  = 


AL 


(85) 


Equations  (68),  (70),  (73),  and  (74)  are  then  substituted  into  (75)  to  obtain: 
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(86) 


■'  dJ,(s)  +  i,^K,+i,,T,(,s)T,{,s) 

The  coefficients  for  the  s  terms  are  also  determined  by  the  MATLAB  code  in  Appendix 
E.  Input  impedance  calculations  are  performed  for  both  the  hard  switching  and  soft 
switching  converter  options.  Note,  Equation  (86)  holds  for  the  specific  multi-loop 
control  with  feedforward  vmder  consideration. 

2.  Determination  of  Filter  Outpnt  Impedance 

Figure  IV-3(a)  is  a  representation  of  the  simplified  network  looking  back  through 
the  filter  to  the  source.  The  network  which  is  used  to  find  the  Thevenin  equivalent 
impedance  is  shown  in  Figure  IV-3(b).  By  appropriately  combining  series  and  parallel 
impedances: 


(a)  (b) 

Figure  rV-3  Filter  Output  Impedance  (a)  Simplified  System  Diagram,  (b)  Thevenin 

Equivalent  Impedance  Network 
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3.  Impedance-Based  Stability  Criteria 

The  type  of  stability  used  for  this  particular  analysis  is  Bounded  Input  Bounded 
Output  (BIBO)  stability.  Per  reference  [23],  “A  system  is  stable  if  every  bounded  input 
excites  a  bounded  output.  Otherwise  the  system  is  said  to  be  imstable.”  A  system  is 
determined  to  be  stable  if  no  poles  of  the  open-loop  transfer  function  are  located  in  the 
open  right-hand  s-plane.  Referring  back  to  Equation  76,  it  is  known  that  is  designed 
to  be  stable  through  the  placement  of  poles  far  in  the  left-hand  s-plane  with  the  effective 
PID  (multi-loop)  controller.  The  filter  transfer  ftmction  (/z^ )  is  also  always  BIBO  stable 

(it  is  a  strictly  passive  network);  therefore,  the  product  of  these  two  terms  will  conversely 
be  BIBO  stable.  As  a  result  of  this,  stability  of  the  overall  system  can  be  determined  by 
analyzing  the  denominator  of  Equation  76  using  the  Nyquist  stability  theorem.  In  order 
for  the  system  to  be  stable,  the  Nyquist  stability  theorem  requires  that  the  Nyquist  plot  of 
the  system  encircle  the  (-1  +  Oj)  point  counter-clockwise  the  number  of  times  equal  to  the 
number  of  poles  in  the  open  right-hand  s-plane.  Since  there  are  no  system  poles  in  the 
open  right-hand  s-plane  in  this  case,  no  encirclements  of  the  (-1  -i-Oj)  point  should  occur. 
To  further  ensure  that  stability  exists,  the  Nyquist  plot  of  the  system  is  confined  to  the 
unit  circle: 

<1  for  -CO  <0X0:^  (88) 

With  these  requirements  in  place,  BIBO  stability  is  guaranteed. 

4.  Evaluation  of  the  Analyzed  System 

A  MATLAB  file  was  written  to  examine  the  stability  of  the  system  using  various 
filter  parameters.  This  code  is  included  in  Appendix  E.  The  Nyquist  plot  using  the  hard 
switching  model  together  with  filter  parameters  w//and  =8265  //F, 

and  a  source  resistance  estimated  at  =  0.01  Q  is  shown  in  Figure  IV-4.  It  can  be 


36 


seen  that  the  Nyquist  contour  is  positioned  within  the  unit  circle.  An  examination  of  the 
poles  shows  that  there  are  no  open  right-hand  s-plane  poles.  It  is  then  confirmed  to  be 
BIBO  stable. 

NyquM  Olftomm* 


Figure  IV-4  Nyquist  Plot  of  BIBO  Stable  System 


The  result  of  lowering  the  filter  capacitor  to  6000  /jF  is  shown  in  Figure  IV-5. 
The  Nyquist  contour  encircles  (-1  +0j)  twice  now.  There  are  still  no  right-hand  s-plane 
poles;  therefore,  BIBO  stability  is  not  assured. 

A  noteworthy  anomoly  was  encoimtered  as  a  result  of  the  method  used  to 
calculate  the  controlled  converter’s  poles  and  zeros.  There  exists  a  pole-zero  cancellation 
that  initially  seemed  not  to  occur.  In  fact,  it  appeared  that  the  system  had  a  RHP  pole. 
Further  investigation  revealed  that  the  cause  was  roundoff  errors  induced  over  multiple 
intermediate  steps  in  the  calculations. 

To  this  point,  two  small-signal  models  have  been  evaluated  to  establish  stability 
criteria  for  the  dc  link  filter.  In  the  next  chapter,  Lyapunov-based  large-signal  stability 
contraints  are  applied  to  a  basic  system  model. 
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Ima^Aicis 


MyQu^t  Diagrioma 


Figure  IV-5  Nyquist  Plot  of  Non-BIBO  Stable  System 
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V.  LYAPUNOV-BASED  LARGE-SIGNAL  STABILITY  CONSTRAINTS 


The  previous  stability  criterion,  specifically  the  small-signal  models,  apply  to 
small  disturbances  in  the  system.  In  Reference  [21],  the  authors  examine  a  stricter  set  of 
criteria  which  act  to  bound  the  source  resistance  for  a  given  set  of  filter,  source  voltage, 
and  load  power  parameters.  They  note  that  when  small-signal  models  are  used,  the 
stability  margins  that  are  used  for  the  final  design  are  forced  to  be  very  conservative  to 
ensure  stability  during  a  large  disturbance.  The  two  criteria  specifically  examined  are  the 
equilibrium  point  criterion  and  the  mixed  potential  criterion.  This  thesis  then  applies 
these  criteria  to  representative  values  being  considered  for  the  DC  ZEDS  project. 


A.  EQUILIBRIUM  POINT  CRITERION 


The  system  examined  is  the  same  as  in  Figure  III-l .  Equations  (43)  and  (44)  are 
the  equations  that  describe  the  states  of  the  system.  The  equilibrium  points  of  these 
equations  are  found  by  setting  the  derivatives  of  and  with  respect  to  time  to  zero: 


and 


(89) 


(90) 


where  capital  letters  are  used  to  distinguish  steady-state  quantities  from  instantaneous 
values. 


In  order  for  the  above  equilibrium  points  to  be  stable,  they  must  be  positive  and 


real.  For  the  operating  point  to  be  real,  must  be  smaller  than  the 
This  results  in  the  condition: 


term. 


V  ^  y 
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es 


(91) 


K  < — r 
'  4P 


The  autiiors  in  [21]  argue  that  since  this  is  a  hard  requirement,  it  is  prudent  to  apply  a 
safety  factor  of  two  when  using  this  criteria.  The  combined  source-filter  resistance  must 
satisfy: 


R 


e,warg 


(92) 


B.  MIXED  POTENTIAL  CRITERION  (SINGLE  CONSTANT  POWER 
LOAD) 


The  authors  next  examine  the  mixed  potential  criterion,  which  was  originally 
developed  specifically  for  non-linear  systems.  The  following  mixed  potential  equation  is 
found  using  the  same  system  as  that  used  for  the  eqmlibriiun  point  criterion: 

U  =  —Reil  +  +  irises  "  V*,  )  (93) 

^<fc/ 

Lyapunov-based  stability  conditions  are  then  placed  on  this  system  which  results  in  the 
following  requirement  to  ensure  that  the  system  response  to  a  large  disturbance  is  not 
oscillatory  [see  21]: 


Since  this  is  not  a  hard  value,  no  safety  margin  is  required  for  this  criteria.  Thus,  the 
design  specification  on  the  source  plus  filter  resistance  combines  Equations  (92)  and  (94) 
to  yield: 


(95) 


Thus,  a  low  equivalent  source  resistance  is  attractive  from  the  standpoint  of 
regulation  and  efficiency,  but  it  in  turn  may  require  a  large  to  ensure  system  stability. 
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The  limitation  of  this  analysis  is  that  a  fairly  simple  representation  of  the  converter  is  still 
employed.  In  particular,  little  can  be  deduced  as  to  the  impact  of  control  bandwidth  on 
the  stability  conditions.  Also,  it  relies  on  a  very  basic  representation  of  the  synchronous 
machine  and  rectifier.  Since  the  criteria  does  not  guarantee  a  stability  margin,  simulation 
must  be  used  to  assess  transient  response  characteristics.  The  attendent  simplicity, 
however,  does  provide  a  convenient,  simple  tool  which  can  be  used  to  provide  baseline 
conclusions. 


C.  MIXED  POTENTIAL  CRITERION  (MULTIPLE  CONSTANT  POWER 
LOADS) 


In  the  dc  power  distribution  system  envisioned  for  DD-21 ,  there  are  multiple 
converters  and  loads,  each  with  their  own  link  filters,  operating  in  parallel  supplied  by  the 
same  source.  This  may  represent  a  generator-rectifier  supplying  a  dc  bus  which,  in  turn, 
delivers  power  to  the  various  zones  within  the  ship  through  separate  regulated  dc-dc 
converters  with  associated  input  filters.  It  is  important  to  understand  how  the  additional 
converters  affect  the  constraints  placed  on  the  filter  design.  The  system  representation  of 
a  multiple  constant  power  load  system  is  shown  in  Figure  V-1.  The  operating  point 
criterion  is  applied  to  this  configuration  to  yield: 


J/2 

^ _ __es _ 

'  4(P1  +  P2  +  P3) 


(96) 


For  simplicity,  the  case  is  considered  where  each  of  the  filters  and  constant  power 
loads  are  identical,  therefore  =  L-^  =  L ,  =  C2  =  =  C ,  and 

Pl  =  P2  =  P3  =  P .  Then  the  Lyapimov  stability  criterion  requires  establishing  bounds 
on  the  norm  of  [see  21]: 


K  = 


-4l 


-2  2 
1  -2 
0  1 


(97) 
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If  an  upperbound  of  the  norm  of  K  is  constrained  to  be  less  than  unity  for 
stability,  the  following  relationship  for  source  resistance  holds: 

Rg  -n 

Vc 

where  n  is  the  number  of  converters.  For  the  case  of  three  converters  illustrated  in 
Figure  V-1: 

^,>4vi5 


Figure  V-1  Multiple  Constant  Power  Load  System  Representation 


(100) 


The  ACSL  code  of  a  single  load  system  is  contained  in  Appendix  F.  Several 
simulations  of  the  system  were  conducted  to  verify  the  validity  of  the  mixed  potential  and 
equilibrium  point  criteria.  The  values  used  for  the  simulation  were  P*  =  200  kW,  = 

1  lOOV,  =  0.5  mH,  and  Q  =  8265  //  F.  These  parameters  result  in  boimdaries  for 
R  between  0.2460  Q  and  1.5125  Q.  When  an  R  of  0.017  Q  is  used,  which  is  below  the 
lower  limit  defined  by  the  mixed  potential  criterion,  it  can  be  seen  in  Figure  V-2  and 
Figure  V-3  that  the  system  is  stable  for  a  small  5%  disturbance  in  the  system  voltage  but 
is  unstable  for  a  large  25%  disturbance  in  the  system  voltage.  In  both  figures,  the 
disturbance  is  introduced  at  time  0.5  seconds.  The  initial  transients  are  due  to  the  system 
reaching  steady-state  from  the  initial  simulation  conditions,  which  were  approximated  to 
be  close  to  the  anticipated  steady-state  values. 


T  t.  i  mie 


Figure  V-2  System  Response  With  R=0.017  Ohms,  5%  Drop 
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Figure  V-3  System  Response  With  R=0.017  Ohms,  25%  Drop 


Likewise,  when  a  value  for  R  of  0.26  Q  is  used,  it  can  be  seen  in  Figure  V-4  and 
Figure  V-5  that  the  system  is  stable  for  both  small  (5%)  and  large  (25%)  disturbances  in 
the  system  voltage.  Like  the  previous  simulations,  in  both  figures,  the  disturbance  is 
introduced  at  time  0.5  seconds.  The  initial  transients  are  again  due  to  the  system  reaching 
steady-state  from  the  initial  simulation  conditions,  which  were  approximated  to  be  close 
to  the  anticipated  steady-state  values. 


Figure  V-4  System  Response  With  R=0.26  Ohms,  5%  Drop 
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Load  1  Current  Load  1  Voltaoe 
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VI.  INTEGRATED  SYSTEM  REPRESENTATION 


A.  USEFULNESS  OF  AN  INTEGRATED  SYSTEM 

The  models  used  to  analyze  the  stability  criteria  were  all  limited  in  scope  to 
varying  degrees  due  to  necessary  simplifying  assumptions.  A  more  detailed  model  of  the 
system  under  consideration  would  be  useful  to  verify  the  validity  of  stability  criteria 
established  by  different  methods.  Models  of  both  the  synchronous  generator-phase 
controlled  rectifier  pair  and  the  filter-converter  pair  have  been  presented  earlier  in  this 
thesis.  While  useful,  these  models  are  also  limited  in  that  the  dynamics  of  the  load  are 
not  experienced  in  the  generator-rectifier  model  and  the  dynamics  of  the  source  are  not 
experienced  in  the  filter-converter  model. 

B.  INTEGRATED  SYSTEM  MODEL 

A  circuit  diagram  of  the  system  modeled  by  the  ACSL  file  contained  in  Appendix 
G  is  shown  in  Figure  VI- 1 .  A  block  diagram  of  the  model  which  clearly  shows  the 
various  subsystems  represented  in  the  model  and  the  interaction  of  the  system  parameters 
between  the  subsystems  is  shown  in  Figure  VI-2.  The  four  main  subsystems  represented 
in  the  model  are  the  synchronous  machine  and  rectifer,  the  generator  field  regulator,  the 
link  filter  and  the  dc-dc  converter.  The  synchronous  generator/rectifier  section  is  based 
on  the  equations  derived  in  Chapter  II.  With  a  dynamic  system  model,  the  need  for 
voltage  control  of  the  generator  becomes  evident.  The  derivations  for  the  equations 
which  govern  the  synchronous  machine  voltage  regulation  are  contained  in  Appendix  B. 
The  link  filter  and  dc-dc  converter  section  is  identical  to  that  discussed  in  Chapter  IV.  Of 
interest  is  the  fact  that  the  synchronous  generator/rectifier/regulator  portions  are  modeled 
in  per  unit  parameters,  while  the  filter/converter  portions  are  modeled  in  actual 
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component  parameters.  This  is  accoimted  for  in  the  code  by  utilizing  conversion  factors 
on  the  system  variables  passed  between  the  two  sections. 


Idc 


Figure  VI-1  Circuit  Diagram  of  Integrated  System 


Figure  VI-2  Block  Diagram  of  Integrated  System 


Several  simulations  were  run  using  the  model  to  verify  the  validity  of  both  the 
model  and  the  design  criteria  evaluated  in  this  thesis.  The  results  of  an  up-power 
transient  from  10%  rated  power  to  100%  rated  power  at  time  0.52  seconds  and  the 
subsequent  down-power  transient  back  to  10%  rated  power  at  time  0.57  seconds  are 
shown  in  Figure  VI-3.  The  system’s  component  values  were  all  such  that  the  various 
stability  requirements  were  met.  It  is  shown  that  the  system  response  is  stable  and  has  an 
acceptable  response  time. 


Figure  VI-3  Up-Power,  Down-Power  Transient,  Cf  =  8265  uF 
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The  same  up-power,  down-power  transient  with  the  filter  capacitor  value  changed 
to  400  //F  is  shown  in  Figure  VI-4.  While  this  value  of  capacitance  is  enough  to 
maintain  stability  at  the  10%  rated  power  level,  it  is  not  sufficient  to  maintain  stability  at 
100%  rated  power.  This  is  confirmed  by  the  stability  criterion  presented  in  earlier 
chapters.  Large,  sustained  oscillations  are  obtained  on  the  dc  bus.  It  is  important  to  note 
that  the  dc-dc  converter  controller  maintains  the  load  voltage  at  the  desired  reference 
voltage  even  with  the  sustained  oscillations  present  at  the  converter  input. 


CO 


Figure  VI-4  Up-Power,  Down-Power  Transient,  Cf  =  400  uF 
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A  couple  of  nuances  that  were  encountered  while  conducting  these  system  studies 
deserve  further  investigation.  First,  optimization  of  the  numerical  integration  time  step 
should  be  performed  to  uncover  at  what  point  the  system  becomes  numerically  unstable 
and  at  what  point  it  becomes  inaccurate.  Further,  the  steady-state  duty  cycle  did  not 
appear  to  settle  at  the  anticipated  value  of  0.88.  It  seemed  to  hone  in  on  a  value 
approximately  5%  higher.  This  anomoly  should  be  explained.  Third,  the  filter  current 
should  be  bounded  so  that  there  are  no  occasions  when  it  can  go  negative.  Finally,  the 
voltage/current  scaling  procedure  should  be  replaced  by  a  more  methodical  scheme. 

The  ACSL  model  of  the  representative  system  has  been  shown  to  be  a  useful  tool 
to  verify  the  results  of  stability  analyses  conducted  by  any  method.  The  next  chapter  will 
present  a  summary  of  the  findings  of  this  thesis  and  discuss  future  work  that  can  be 
performed  to  futher  knowledge  about  this  issue. 
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VII.  CONCLUSIONS 


A.  SUMMARY  OF  FINDINGS 

Two  small-signal  models  and  a  large-signal  model  of  a  representative  dc  electric 
power  distribution  system  are  examined  to  determine  the  design  criteria  for  the  dc  link 
filter  in  the  system.  The  Lyapunov-based  large-signal  stability  criterion  is  fovmd  to  be  the 
most  conservative  of  the  three  analyzed.  The  next  most  conservative  is  the  simplified  dc 
link  model  with  the  Hurwitz  criteria  for  stability  applied.  The  least  conservative  is  the 
application  of  the  Nyquist  criteria  to  the  impedance-based  model.  Simulations  of  the 
different  simplified  models  were  used  to  illustrate  the  validity  of  the  filter  design 
limitations  defined  by  each  criterion. 

When  component  values  that  meet  the  simplified  model  stability  criteria 
examined  are  used,  detailed  model  simxilations  show  that  these  criteria  are  conservative 
enough  to  be  considered  valid  for  hardward  design.  A  supplemental  finding  worthy  of 
comment  is  the  region  of  stability  outside  of  the  criteria  bounds  established  upon 
investigation  of  the  detailed  filter-converter  model.  The  region  exists  as  a  result  of  the 
feedback  controller  utilized  by  the  converter  to  control  the  duty  cycle. 

B.  FUTURE  WORK 

This  thesis  motivates  the  investigation  of  other  more  detailed  aspects  of  stability 
concerns  in  a  DC  ZEDS  design.  First,  further  investigation  into  the  stability  concerns 
involved  with  operating  multiple  controlled  converters  in  parallel  can  be  undertaken.  The 
simplified  models  used  in  this  thesis  ignored  the  effects  of  controlled  converters.  The 
relationship  between  control  bandwidth  and  stability  should  be  investigated.  An 
additional  issue  is  raised  when  the  load  on  the  dc  system  is  an  inverter  vice  a  converter. 
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The  input  characteristics  of  a  tightly-regulated  inverter  differ  from  those  of  a  dc-dc 
converter. 

This  thesis  has  focussed  on  meeting  the  system  stability  requirements  by  proper 
selection  of  the  passive  components  in  the  device  input  filters.  Further  study  should  be 
performed  with  regards  to  introducing  nonlinear  compensators  into  the  system.  These 
compensators  would  preclude  the  need  for  large  capacitor  banks  which  are  both  difficult 
to  troubleshoot  and  have  marginal  reliability.  Finally,  the  active  compensation 
techniques  may  be  applied  to  the  synchronous  machine  field  regulator,  the  bus  voltage 
phase-controlled  rectifier,  or  with  regards  to  the  dc7dc  converter  controls.  Many  research 
opportunities  are  currently  available  in  this  area  and  should  be  explored. 

Factors  that  were  assumed  negligible  in  this  thesis  that  may  need  to  be  examined 
in  the  future  include  the  effects  of  line  inductance  and  parasitic  capacitance  in  a  full  size 
system.  Additionally,  the  effects  of  electro-magnetic  interference  (EMI)  due  to  the  high 
switching  frequencies  used  in  the  converters  was  not  considered. 

In  conclusion,  various  stability  criteria  have  been  examined  to  address  the 
stability  issues  involved  with  having  link  filters  as  an  integral  part  of  a  dc  electrical 
distribution  system.  The  analysis  and  tools  crafted  in  this  thesis  form  a  strong  foundation 
for  establishing  design  specifications  for  the  L-C  filters  in  the  DC  ZEDS  program. 
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APPENDIX  A:  SYNCHRONOUS  MACHINE  -  RECTIFIER  ACSL  MODEL 


A.  .CSL  FILE 


PROGRAM  ! " tues . csl " 


l”Greg  Greseth,  LT,  USN  -  8  Jan  99” 

!  "Last  Revised  -  15  Feb  99 

!  "This  ACSL  file  simulates  a  reduced  order  model  for  a 


synchronous” 

are” 


! "machine  and  three  phase  rectifier. 
! "defined  in  *per  unit*  values." 


The  machine  parameters 


INITIAL  ! "Section  in  which  constants  and  program 


parameters 

are  defined" 

ALGORITHM  ialg  =  5 

! "Runge-Kutta  4th  order 

CONSTANT  tstop  =0.5 

! "stop  time  for 

integration 

n 

MAXTERVAL  maxt  =  l.e-5 

! "max  integration  step 

size” 

MINTERVAL  mint  =  l.e-6 

! "min  integration  step 

size" 

CINTERVAL  cint  =  5.e-4 

! "data  communication 

interval" 

CONSTANT  wb  =  377. 

CONSTANT  rkgpri  =  0.0613 
CONSTANT  rkdpri  =0.24 

CONSTANT  rfdpri  =  0.00111 
CONSTANT  J  =  1. 

CONSTANT  T1  =  1. 

CONSTANT  B  =  0. 

CONSTANT  P  =  2. 

CONSTANT  Xls  =0.08 

CONSTANT  Xmq  =  1.00 

CONSTANT  Xlkqpri  =  0.330 

CONSTANT  Xmd  =1.77 

CONSTANT  Xlfdpri  =  0.137 

CONSTANT  Xlkdpri  =  0.334 

CONSTANT  rs  =  0.00515 

! "number  of  poles" 

CONSTANT  exfd  =  1. 

! "CONSTANT  beta  =  0.523599" 

! "CONSTANT  alpha  =0." 

CONSTANT  alfahat  =0.0 

!"2.15" 

CONSTANT  del  =  0.0 

CONSTANT  PI  =  3.141593 

CONSTANT  radl20  =  2.094395 
CONSTANT  rad60  =  1.047198 
CONSTANT  rad30  =  0.523599 

1 "0.9" 
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u  =  0.05  ["initialize  u  for  use  in  first  findu  function 

call" 

’"define  the  subtransient  reactances" 

Xqpripri  =  Xls+ (Xmq*Xlkqpri) / (Xlkqpri+Xmq) 

Xdpripri  = 

Xls+ (Xmd*Xlfdpri*Xlkdpri) / (Xmd*Xlfdpri+Xmd*Xlkdpri& 

+Xlkdpri*Xlfdpri) 

Vfdpri  =  (exfd*rfdpri) /Xmd 

COMMON  /  forsub  / 

rad30, rad60, radl20^Xqpripri, Xdpripri, del, PI 
END  1 "of  initial  section" 

DYNAMIC 

TERM? (t .ge. tstop,  ’Termination  on  time  limit’) 

DERIVATIVE 

["state  equations  for  psikqpri,  psikdpri,  psifdpri,  wr 

and  thetar" 

psikqpri  =  INTEG (-wb*rkqpri*ikqpri ,  0.0)  ["0.049" 

psikdpri  =  INTEG (~wb*rkdpri*ikdpri,  1,0)  ! "0.293" 

psifdpri  =  INTEG (-wb^rfdpri^ifdpri+wb^Vfdpri,  1.077) 

["0.460" 

[ "wrm  =  INTEG( ( (l./J) * (Te-Tl-B*wrm) ) ,  0.)" 

[ "relationship  between  wrm  and  state  wr 

[ "wr  =  (P/2.)*wrm" 
wr  =  377. 

thetar  =  INTEG(wr,  0.) 

["thetar  =  mod {theta, rad60 )  +  radGOd" 

[ "electrical  torque  equation" 

I  iripg  _  ^  ^  /pj  /  (4  ^  *wb)  )  *  (psidsr*iqsr-psiqsr*idsr )  " 

Te  =  psidsr*iqsr-psiqsr*idsr 

["solve  first  for  ikdpri,  ifdpri  and  ikqpri" 

ikdpri  =  (1 .  /  (  (Xlkdpri-fXmd)  *  (Xlfdpri+Xmd) - 
Xmd**2 .  )  )  *  (  (Xlfdpri+Xmd)  & 

* (psikdpri-Xmd*idsr) -Xmd* (psifdpri-Xmd^idsr ) ) 


56 


ifdpri  -  (l./(  {Xlkdpri+Xmd)  *  (Xlfdpri+Xmd)-Xitid**2. )  )  *  ( 

Xmd* (psikdpri& 

-Xmd*idsr) + (Xlkdpri+Xmd) * (psifdpri-Xmd^idsr ) ) 

ikqpri  =  psikqpri/ (Xlkqpri+Xmq) - 
{Xmq*iqsr) / (Xlkqpri+Xmq) 

! "equations  to  eliminate  psiqsr  and  psidsr" 

psiqsr  =  Xqpripri*iqsr+psiqpripri 

psidsr  =  Xdpripri^idsr+psidpripri 

! "define  the  subtransient  flux  linkages" 

psiqpripri  =  (Xmq/ (Xmq+Xlkqpri) ) *ps ikqpri 

psidpripri  = 

( (Xmd*psikdpri) /Xlkdprit (Xmd*psifdpri) /Xlfdpri) & 

/ {1 , +Xmd/Xlfdpri+Xmd/Xlkdpri) 

! "define  Eqpripri  and  Edpripri" 

Eqpripri  =  (wr*psidpripri ) /wb 

Edpripri  =  - (wr*psiqpripri ) /wb 

{"define  phi  to  use  to  determine  alpha" 

phi  =  ATAN2 (Edpripri, Eqpripri) 

1 "call  function  to  determine  beta" 

beta  =  findbeta (psidpripri, psiqpripri, alfahat) 

! "define  alpha" 

alpha  -  beta  -  phi 

! "desired  reduced  order  voltage-behind-reactance 

equations" 

Vqsr  =  rs*iqsr+ (wr/wb) *Xdpripri*idsr+Eqpripri 

Vdsr  =  rs*idsr- (wr/wb) *Xqpripri*iqsr+Edpripri 

! "equation  for  average  DC  voltage  out  of  rectifier" 

Vdc  =  (3.*sqrt (3. ) /PI)* (wr/wb) * (sqrt (psiqpripri**2& 
+psidpripri**2) ) * (cos (alpha) ) - 
(3. /PI) * (wr/wb) * ( ( (Xqpripri+Xdpripri) /2. ) *dcl& 

+ (Xdpripri-Xqpripri) *dcl*sin (2 . *beta+rad30) ) 

{"need  to  define  K1  for  use  in  subroutines" 

Kone  =  sqrt (3 . ) * (- 

psiqpripri*sin (beta) +psidpripri*cos (beta) ) & 
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+  {  (Xdpripri-Xqpripri)  *cos  (2 .  *beta+racil20 )  -0 . 5*  (Xqpripri+Xdpripri)  )  *dcl 

! "call  function  to  determine  u  using  ’half-it’  method" 
u  =  findu (Kone,psiqpripri,psidpripri) 

! "now  with  the  new  u,  we  can  determine  conduction 

currents" 


-cos (beta+u+rad60 

-sin (beta+u+rad60 

method" 


iqsrcond  =  ( (2 . *sqrt {3, ) ) /PI) *dcl* (cos (beta+radl20) & 

)  ) 

idsrcond  =  ( (2 . *sqrt (3 . ) ) /PI ) *dcl* (sin (beta+radl20) & 

)) 

! "calculate  the  commutating  currents  using  Simpson 
! "first,  define  the  5  different  values  of  thetar  to  be 


used" 


thetarl  =  beta+rad60 
thetar2  =  beta+rad60+u/4 . 
thetar3  =  beta+rad60+u/2 • 
thetar4  =  beta+rad60+ (3 . *u) /4 . 
thetarS  =  beta+rad60+u 

’"call  each  q  and  d  subroutine  5  times  to  calculate 

components" 

PROCEDURAL  (Q1  =  thetarl, Kone, psiqpripri, psidpripri) 
CALL 

simpsonq (thetarl , Kone, psiqpripri, psidpripri, Q1 ) 

END  ! "of  procedural" 

PROCEDURAL  (D1  =  thetarl, Kone, psiqpripri, psidpripri ) 
CALL 

simpsond (thetarl , Kone, psiqpripri, psidpripri,  D1 ) 

END  ! "of  procedural" 

PROCEDURAL  (Q2  =  thetar2 , Kone, psiqpripri, psidpripri ) 
CALL 

simpsonq ( thetar2 , Kone, psiqpripri , psidpripri,  Q2 ) 

END  ! "of  procedural" 

PROCEDURAL  (D2  =  thetar2, Kone, psiqpripri, psidpripri) 
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CALL 

siinpsonci(thetar2,Kone,psiqpripri,psidpripri,  D2) 

END  1 ”of  procedural” 

PROCEDURAL  (Q3  =  thetarS, Kone, psiqpripri, psidpripri) 
CALL 

simpsonq (thetarS, Kone, psiqpripri, psidpripri,  Q3) 

END  ! "of  procedural” 

PROCEDURAL  (D3  =  thetar3, Kone, psiqpripri, psidpripri) 
CALL 

simpsond (thetar3, Kone^psiqpripri^psidpripri,  D3) 

END  ! "of  procedural" 

PROCEDURAL  (Q4  =  thetar4 , Kone, psiqpripri, psidpripri ) 
CALL 

simpsonq  (thetar  4 ,  Kone,  psiqpripri,  psidpripri,  Q4 ) 

END  ! "of  procedural” 

PROCEDURAL  (D4  =  thetar4 , Kone, psiqpripri^ psidpripri ) 
CALL 

simpsond  (thetar4 ,  Kone,  psiqpripri,  psidpripri,  D4 ) 

END  !"of  procedural” 

PROCEDURAL  (Q5  =  thetarS, Kone, psiqpripri, psidpripri) 
CALL 

simpsonq (thetarS, Kone,psiqpripri,psidpripri, Q5) 

END  i "of  procedural” 

PROCEDURAL  (D5  =  thetarS , Kone, psiqpripri, psidpripri) 


CALL 

simpsond (thetarS,  Kone,  psiqpripri^psidpripri,  D5) 

END  !"of  procedural” 

! "now  add  all  the  components  together" 
iqsrcom  =  (u/ ( PI*4 . ) ) * (Ql+4 . *Q2+2 . *Q3+4 . *Q4+Q5 ) 
idsrcom  =  {u/ (PI*4 . ) ) * {Dl+4 . *D2+2 . *D3+4 . *D4+D5) 


average" 


1 "add  the  commutating  and  conduction  components  to  get 
!"iqsr  and  idsr” 
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iqsr  =  iqsrcom+iqsrcond 


idsr  =  idsrcom+idsrcond 


END  ! "of  derivative" 


END  ! "of  dynamic" 

END  ! "of  program" 

SUBROUTINE  simpsonq (thetarsimp, Kone, psiqpripri, psidpripri, iqsrout ) 


*  subroutine  to  calculate  components  for  simpson  method  of  determining 

*  the  q  commutating  currents 

* 


COMMON  /  forsub  /  radSO, rad60, radl20, Xqpripri, Xdpripri, del, PI 

* 

*  declare  variables 

* 

real  ias,  parti,  part2,  part3,  Kone,  iqsrout,  del 

* 

*  break  the  equation  for  ias  down  into  manageable  parts 

partl=Kone-sqrt (3. ) * (psiqpripri*cos (thetarsimp+radSO) +psidpripri 
+  *sin (thetarsimp+radSO ) ) 

★ 

part2= (Xqpripri-Xdpripri) *cos (2 . *thetarsimp+rad60) +Xqpripri 
+  +Xdpripri 

* 

part3=( (Xqpripri-Xdpripri)*cos(2.*thetarsimp+radl20)+0.5 
+  * (Xqpripri+Xdpripri ) ) *dcl 

★ 

*  combine  the  3  parts  to  determine  ias 

* 

ias= (partl/part2) - (part3/part2) 

* 

*  now  solve  for  iqsr  at  the  inputed  thetar 

•A- 

iqsrout=( (2.*sqrt (3. ) ) /3. ) * (ias*cos {thetarsimp+rad30) -del 
+  *sin (thetarsimp) ) 

■k 

end 

SUBROUTINE  simpsond (thetarsimp, Kone, psiqpripri, psidpripri, idsrout) 

★ 

*  subroutine  to  calculate  components  for  simpson  method  of  determining 

*  the  d  commutating  currents 


COMMON  /  forsub  /  rad30, radGQf radl20, Xqpripri, Xdpripri, del, PI 

*  declare  variables 

* 

real  ias,  parti,  part2,  part3,  Kone,  idsrout,  del 
*  break  the  equation  for  ias  down  into  manageable  parts 
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partl=Kone-sqrt (3. ) * (psiqpripri*cos (thetarsimp+rad30) +psidpripri 
+  *sin {thetarsimp+rad30) ) 

* 

part2={Xqpripri-Xdpripri) *cos (2.*thetarsimp+rad60) +Xqpripri 
+  +Xdpripri 

* 

part3= ( (Xqpripri-Xdpripri) *cos (2 . *thetarsimp+radl20) +0 . 5 
+  * (Xqpripri+Xdpripri) ) *dcl 

★ 

*  combine  the  3  parts  to  determine  ias 

* 

ias= (partl/part2) - (part3/part2) 

*  now  solve  for  idsr  at  the  inputed  thetar 

★ 

idsrout=  (  (2  .  *sqrt  (3. )  )  /  3 . )  *  (ias* sin  (thetarsimp4-rad30 )  +dcl 
+  *cos (thetarsimp) ) 

* 

end 

REAL  FUNCTION  findu (Kone, psiqpripri, psidpripri ) 


*  function  subprogram  to  calculate  u  using  the  ’half-it*  method 

★ 

COMMON  /  forsub  /  rad30, rad60, radl20, Xqpripri, Xdpripri, del,  PI 

* 

*  declare  variables  to  be  used 

* 

logical  condition 

real  utemp,  fcom,  fpricom,  Kone,  del,  ul,  u2,  u 

* 

ul=0.0 

u2=rad60 

u=0.0 

f com=Kone-sqrt (3 . ) * (“psiqpripri*sin (beta+u) +psidpripri 
+  *cos (beta+u) ) - ( (Xqpripri-Xdpripri) *cos (2 . *beta+2 . *u 

+  ~radl20) +0.5* (Xqpripri+Xdpripri) ) *dcl 

10  if  (ABS(ul-u2)  .gt.  0,001)  then 
u  =  0.5*(ul+u2) 

fcom=Kone-sqrt (3. ) * (-psiqpripri*sin (beta+u) +psidpripri 
+  *cos (beta+u) ) - ( (Xqpripri-Xdpripri) *cos (2 . *beta+2 . *u 

+  -radl20) +0 . 5* (Xqpripri+Xdpripri) ) *dcl 

if  (fcom  .le.  0.00)  then 
ul=u 
else 
u2=u 
endif 
goto  10 
endif 
findu=u 
END 

REAL  FUNCTION  findbeta (psidpripri, psiqpripri,  alfahat) 


*  function  subprogram  to  calculate  beta 
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COMMON  /  forsub  /  radSO, rad60, radl20,Xqpripri,Xdpripri^  del, PI 


*  declare  variables  to  be  used 

* 

real  betatop,  betabot,  betamin,  alfahat,  fbmin 

* 

betatop  =  PI 
betabot  =  -PI 
betamin  =  0.0 

* 

*  calculate  the  function  to  determine  if  beta  is  greater  than  betamin 

* 

fbmin  =  SQRT (3 . ) * (psiqpripri^cos (betamin) +psidpripri*sin (betamin) ) 
+  +2 . *dcl* (Xqpripri-Xdpripri) *sin (2 . *betamin-rad60) 

10  if  (ABS (betatop-betabot)  .gt.  0.001)  then 
betamin  =  (betatop+betabot ) /2 . 

fbmin  -  SQRT (3 .)* (psiqpripri*cos (betamin) +psidpripri 
+  *sin  (betamin) ) +2 . *dcl* (Xqpripri-Xdpripri) *sin (2 . ^betamin 

+  -rad60) 

if  (fbmin  .le.  0.0)  then 
betabot  =  betamin 
else 

betatop  =  betamin 
endif 
goto  10 
endif 

★ 

*  now  determine  beta  based  on  the  above  determination 

* 

findbeta  =  alfahat  +  betamin 
END 


B.  .CMD  FILE 


1 "tues . cmd" 

!"Greg  Greseth  -  15  Feb  1999" 

! "This  file  plot  the  output  of  a  synchronous  generator,  phase- 
controlled" 

1 "rectifier  model." 


s  strplt  =  .t. 
s  calplt  =  .f. 

s  devplt  =  1 

s  ppoplt  =  .f. 
s  xinspl  =  6 
s  weditg  -  .f. 


!  "one  variable  per  x-axis" 

!  "1  for  X- windows" 

!  "5  for  postscript" 

!  "true  rotates  plot  90  deg" 

!  "x-axis  plot  units" 

!  "false  suppresses  data  write" 


s  wesitg  =  . f . 
s  nrwitg  =  .f. 

s  alcplt  =  .f. 
s  tstop  =  40. 
s  inaxt=1.0e-4 
s  cint=1.0e-3 


!  "each  time  SCHEDULE  occurs" 

!  "true  enables  accumulation  of  data" 
!  "after  a  CONTIN" 

!  "false  makes  all  plots  dark" 


prepare 

t, Vdc, Te,psikqpri,psikdpri,psifdpri,  iqsr, idsr,u,  thetar,beta, alpha, phi 

PROCED  startup 
s  ialg=5 

s  devplt  =  1 
start 

s  pit  =  1 

plot  Vdc,Te,psikqpri,psikdpri 
s  pit  =  2 

plot  psifdpri, iqsr, idsr, u 
s  pit  =  3 

plot  beta,  phi,  alpha 

END 

PROCED  testalfa 

s  tstop  =  60 
s  ialg=5 
start 

d  Vdc, Te,psikqpri,psikdpri, psifdpri, iqsr, idsr, u, beta, phi, alpha 

END 
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APPENDIX  B:  DERIVATION  OF  SYNCHRONOUS  MACHINE  VOLTAGE 

REGULATION 


The  reduced-order  (stator  transients  ignored)  electrical  equations  of  a  3-phase 
synchrounous  machine  with  a  single  q  and  d-axis  damper  winding  may  be  expressed  by: 


(B-1) 

(B-2) 

(B-3) 

^qdr  —Hi  J  i 

(B-4) 

where: 


i 


‘fd 


(B-5) 


(B-6) 


With  the  stator  currents  referenced  into  the  machine  and  assuming  a  purely  resistive  load, 
the  terminal  voltages  may  be  expressed  solely  in  terms  of  the  terminal  currents: 


where: 


~~  ^load  ^ 

0  ~  ^load  . 


(B-7) 


(B-8) 


Upon  substituting  (B-7)  into  (B-1)  and  solving  for  we  get: 

'p;^=-i-'(l-Z)r;^  (b-9) 

Further,  (B-9)  can  be  substitued  into  (B-3)  so  that  the  stator  currents  may  be  solved  for  in 
terms  of  the  rotor  currents: 


65 


(B-10) 


^  qds  ^  ^  qdr 


where: 


l  =  [-b-\a-k)-fYg 

Equation  (B-10)  is  then  plugged  into  (B-4)  to  yield: 

and  solving  for  i  : 

i'qdr=(HL+jr^>'^ 

Finally,  substituting  (B-13)  into  (B-2),  we  get: 


dt 


(B-11) 

(B-12) 

(B-13) 

(B-14) 


and  since  -  0,  the  referred  rotor  voltage  vector  may  be  re-expressed  as: 


=  NV,,  = 


fd 


(B-15) 


The  output  of  the  machine  is  the  terminal  voltages.  Equations  (B-7),  (B-10),  and 
(B-13)  may  be  combined  to  express  the  terminal  voltages  in  terms  of  the  rotor  states: 

V;^=KL(HL+Jy^%,^  (B-16) 

Once  we  assume  that  the  turbine  governor  is  holding  constant,  equations  (B-14)  and 
(B-16)  become  a  set  of  linear  ordinary  differential  equations.  However,  the  peak  value  of 
the  terminal  voltages  is  related  to  and  by: 

(B-n) 


Therefore,  to  consider  a  single-input  (  ),  single-output  (F^^  ^* )  system,  the  output 

relationship  must  be  linearized  about  an  operating  point: 

^Vkpk=^^^qds  (B-18) 

where: 
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(B-19) 


Finally,  the  state-space  model  of  the  generator  is  given  by: 

I (B-20) 

(B-21) 

where: 

I2  =-D''C(^r  +  7)-’  (B-22) 

\  =  D-'N  (B-23) 

C2=MKL(Hl  +  jy'  (B-24) 

A  final  point  lies  in  determining  the  operating  point  values  for  (B-19).  This  can 
be  accomplished  assuming  that  for  rated  load,  Vjy  has  been  selected  to  provide  the  rated 
terminal  voltage.  In  the  steady  state: 


W"  -  -A  "’S  V' 

^  qdro  -^2  ^1*  Jdo 

(B-25) 

Substituting  (B-25)  into  the  output  equation  (B-16),  we  find: 

11 

-KL{HL+jrA{^B,V',, 

(B-26) 

or: 

=  AVfd 

(B-27) 

where  is  a2-by-l  vector: 

7  _r^3(i)' 

—  — 

(B-28) 

_^3  (2)_ 

Thus,  from  (B-17)  and  (B-27): 

(B-29) 

Therefore,  for  a  given  load  resistance  (  A” ),  the  appropriate  field  voltage  to  give  the  rated 
terminal  voltage  is  found  from: 
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1 


(B-30) 


p,(lf+A,(2y 


V 

'  qds,pk,ral 


This  value  may  then  be  substitued  back  into  (B-27)  to  yield  the  desired  operating  point 
values. 


At  this  point,  the  design  of  the  field  regulator  can  commence.  To  illustrate  the 
need  for  compensation,  the  open-loop  eigenvalues  of  A2  fi:om  Equation  B-20  can  be 
calculated  assuming  a  load  resistance  of  1 .0  pu.  These  are  found  to  be: 


=-19.71 
/I2  =  -198.61 


A3  =  -0.54 


The  A3  eigenvalue  introduces  a  fairly  long  recovery  time.  In  addition,  open-loop, 
the  operator  would  need  to  know  the  load  in  order  to  apply  the  correct  to  maintain  the 

desired  .  For  instance,  with  =  1.0,  the  required  field  voltage  to  yield  F^^  = 

1.0  is: 


=  0.0012789 
or 


4=2.039=^=!-^^;,  (B-31) 

V 

As  the  load  power  is  decreased,  this  excitation  voltage  must  be  lowered  to  maintain 
^qds,pk  ~  1  O-  For  the  above  condition,  the  required  stator  voltage  operating  points  are 
found  to  be: 

F;,  =  0.6813  pu 
FI  =  0.7320  pu 

The  system  can  be  linearized  about  these  operating  points  and  a  transfer  fimction  derived 
between  AF^^  and  AF^  : 
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(B-32) 


Gis)  = 


^  219.7(5 +  271)(5  + 27.7) 
AV'j  (5  +  19.7)(j  +  199)(j  +  0.54) 


The  desired  control  loop  is  then  illustrated  in  Figure  B-1.  The  circuit  can  be 
analyzed  to  yield  the  closed-loop  transfer  function: 


and  the  steady-state  error: 


^V^pk  n.n^+d.d^ 


^  d.d^V;lpk 


where  n^,  rig,  d^  and  dg  are  polynomials  in  s. 


(B-33) 


(B-34) 


1 

<1 

d^  Cj  5  + 1 

1 

5 

Figure  B-1  Field  Regulator  Control  Loop 


The  control  gains  Cj  and  must  be  selected  to  yield  satisfactory  steady-state 

error  and  acceptable  closed-loop  pole  locations.  Plugging  in  the  transfer  function 
elements  into  Equation  B-34  yields: 

1  ^ 


E,= 


1  +  799C, 


(B-35) 


If  the  output  voltage  is  to  be  within  0.5%  of  the  commanded  value,  this  imposes  the 
constraint: 

c,  >0.223  (B-36) 
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If  it  is  assumed  that  Cj  =  0.25  for  convenience,  may  be  selected  based  on  the  closed- 
loop  poles  of  Equation  B-33,  Note,  the  zeros  of  Equation  B-33  are  fixed  by  the  zeros  of 
Hg  and  once  c,  is  fixed,  the  steady-state  error  is  fixed.  Table  B-1  illustrates  how  the 
closed-loop  poles  vary  as  is  changed. 


^2 

Pole  1 

Pole  2 

Poles  3  &  4 

1000 

-199 

-19.7 

-0.2697j0.18 

10 

-199 

-19.9 

-0.2357j3.23 

1 

-199 

-21.0 

-0.0847j9.97 

0.1 

-200 

-25.1 

-2.067j28.7 

0.01 

-212 

-30.6 

-38.237j70.23 

0.002 

-457 

-32.6 

-1147j28.3 

0.001 

^  -948 

-33 

-89.1  &  -149 

Table  B-1  Closed-loop  Poles  Versus  Gain  C2 


As  the  gain  is  decreased,  the  poles  3  &  4  move  from  slow,  lightly  damped 
locations  to  positions  comfortably  in  the  left-half  plane.  An  attractive  compensator 
transfer  function  is  given  by: 


or 


AF, 


fit 


0.25 

0.0025 -Hi 


(B-37) 


Ac 


400 


AC,* 


0.0025 -Hi 

One  final  point  to  understand  is  that  this  design  is  valid  for  one  operating 
condition  =1.0  pu).  As  illustrated  by  Table  B-2,  the  steady-state  value  of 


(B-38) 
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must  be  varied  to  maintain  rated  output  voltage.  The  corresponding  transfer  functions 
can  be  evaluated  and  analyzed  for  suitability. 


f 

1 

^xjd  (pu) 

1.0 

2.039 

1.2 

1.771 

1.4 

1.601 

1.6 

1.478 

1.8 

1.388 

2.0 

1.321 

5.0 

1.057 

10.0 

1.015 

Table  B-2  Excitation  Voltage  Required  to  Yield  Rated  Output  Voltage 
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APPENDIX  C:  SIMPLIFIED  DC  LINK  ACSL  MODEL 


A.  NON-LINEAR  .CSL  FILE 


PROGRAM  ! "simplenonlin. csl" 

! "Greg  Greseth  -  8  April  1999" 

! "This  file  contains  the  non-linear  system  equations  for  a  simplified" 
! "DC  link  model.  The  simulated  event  is  a  10%-100%-10%  power 
transient." 


INITIAL 

ALGORITHM  ialg  =  5 
MAXTERVAL  maxt  =  l.e-5 
MINTERVAL  mint  =  l.e-6 
CINTERVAL  cint  =  5.e-3 
CONSTANT  Le  =  0.0005 
CONSTANT  Ce  =  0.008265 
CONSTANT  Re  =  0.01 
CONSTANT  tstop  =5.0 
CONSTANT  P  =  20000. 

CONSTANT  Ves  =  1102. 

END  ! "of  initial 

DYNAMIC 

TERMT(t.ge. tstop,  'Termination  on  time  limit') 
DERIVATIVE 

P_star  =  P  +  180000. *STEP{2.0)  -  180000 . *STEP (3 . 5) 
ir  =  INTEG(-{Vdci/Le)-(Re*ir)/Le+Ves/Le,181.82) 
Vdci  =  INTEG(ir/Ce-{P_star)/(Ce*Vdci),1100.) 

END  ! "of  derivative" 

END  ! "of  dynamic" 

END  !"of  program" 
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B. 


NON-LINEAR  .CMD  FILE 


!  "sirtiplenonlin .  cmd" 

!”Greg  Greseth  -  8  April  1999” 


!  "This  file  plots  the  inductor  current  and  capacitor  voltage” 
!”for  the  non-linear  simplified  DC  link  model” 


s  strplt  =  ,t. 
s  calplt  =  .f. 

s  devplt  =  1 

s  ppoplt  =  .f. 
s  xinspl  =  6 
s  weditg  =  .f. 
s  wesitg  =  .f. 
s  nrwitg  =  .f. 

s  alcplt  =  •f. 


!  ”one  variable  per  x-axis” 


!  ”1  for  X- windows” 

!  ”5  for  postscript" 

!  "true  rotates  plot  90  deg” 

!  "x-axis  plot  units” 

!  "false  suppresses  data  write” 
!  "each  time  SCHEDULE  occurs” 


!  "true  enables  accumulation  of  data” 
!  "after  a  CONTIN” 

!  "false  makes  all  plots  dark” 


prepare  t,ir,Vdci 


PROCED  nonlinear 
start 

plot  ir,Vdci 

END 


C.  LINEAR  .CSL  FILE 


PROGRAM  ! "simplelin. csl" 

!”Greg  Greseth  -  8  April  1999” 

! "This  file  contains  the  linearized  system  equations  for  a  simplified” 
! "DC  link  model.  The  simulated  event  is  a  10%-100%-10%  power 
transient . ” 

INITIAL 

ALGORITHM  ialg  =  5 
MAXTERVAL  maxt  =  l.e-5 
MINTERVAL  mint  =  l.e-6 
CINTERVAL  cint  =  5.e-3 
CONSTANT  Le  =  0.0005 
CONSTANT  Ce  =  0.008265 
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CONSTANT  Re  =  0.01 
CONSTANT  tstop  =1.5 
CONSTANT  P  =  20000. 

CONSTANT  Vdci  =  1100. 

CONSTANT  delta_Ves  =5. 

END  ! "of  initial 

DYNAMIC 

TERMT (t.ge. tstop,  'Termination  on  time  limit') 

DERIVATIVE 

P_star  =  P  +  180000. *STEP(0. 5)  -  180000. *STEP (3 . 0) 

delta_ir  =  INTEG (- {delta_Vdci/Le) - (Re*delta_ir) /Le+delta_Ves/Le, 0) 

delta_Vdci  =  INTEG (delta_ir/Ce+ (P_star*delta_Vdci) / (Ce*Vdci**2) , 5) 

END.  !"of  derivative" 

END  ! "of  dynamic" 

END  ! "of  program" 

D.  LINEAR  .CMD  FILE 


! "simplelin.cmd” 

!”Greg  Greseth  -  8  April  1999" 

!  "This  file  plots  the  linearized  inductor  current  and  capacitor'’ 
'"voltage  for  the  simplified  DC  link  model" 


s  strplt  =  .t. 
s  calplt  =  .f. 

s  devplt  =  1 


s  ppoplt  =  .f. 
s  xinspl  =  6 
s  weditg  =  .f. 
s  wesitg  =  .f. 
s  nrwitg  =  .f. 


!  "one  variable  per  x-axis" 

!  "1  for  X-windows" 

!  "5  for  postscript" 

!  "true  rotates  plot  90  deg” 

!  ”x-axis  plot  units” 

!  "false  suppresses  data  write” 
!  "each  time  SCHEDULE  occurs" 


!  "true  enables  accumulation  of  data" 
!  "after  a  CONTIN" 


s  alcplt  =  .f. 


"false  makes  all  plots  dark" 
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prepare  delta_ir, delta  Vdci 


PROCED  linear 
start 

plot  delta_ir, delta_Vdci 

END 
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APPENDIX  D:  FILTER-CONVERTER  ACSL  MODEL 


A.  .CSL  FILE 


! "hards w. csl” 

! "written  by  John  Ciezki” 

! "modified  by  Greg  Greseth  -  29  April  1999" 

’"last  modified  -  1  May  1999" 

! "This  models  a  DC-DC  converter  (buck  chopper)  with  an  DC  link  filter 
at" 

!"the  input." 


PROGRAM 


INITIAL 


MAXTERVAL  maxt  =  2.0e-7 
MINTERVAL  mint  =  l.Oe-7 


!  "maximum  integration  step  size" 


CINTERVAL  cint  = 
ALGORITHM  ialg  = 

NSTEPS  nstp  = 
CONSTANT  tstop 


l.Oe-3  ’"data  communication  interval" 

5  ! "integration  algorithm" 

!"4“R.K.  2nd,  5— R.K.  '4th" 

1 

=  0.25  !"stop  point  for  integration" 


! "Buck  Parameters" 
CONSTANT  Lbuck  =  l.Oe-3 
CONSTANT  Riser  =0.0 
CONSTANT  Vdidrop  =0.0 
CONSTANT  Vswdrop  =0.0 
CONSTANT  Cbuck  =  2.0e-3 
CONSTANT  Tramp  =  250.0e-6 


"main  inductance" 

"Lbuck  series  resistance" 
"diode  forw  conduction  drop" 
"switch  forward  drop" 

"main  output  capacitance" 
"switching  period  (4kHz)" 


!" - define  and  initialize  the  main  switch  conduction  status" 

LOGICAL  sw 

sw  =  .false.  ! "T  sw  conducts,  F  sw  blocks" 

1 " - define  and  initialize  variable  which  keeps  track  of" 

’"  whether  the  converter  is  in  contin  or  discontin  mode" 
LOGICAL  contmode 
contmode  =  .true. 


’"link  filter  parameters" 
CONSTANT  Cf  =  8265.0e--6 
CONSTANT  Lf  =  0.5e-3 
CONSTANT  Rs  =  0.01 


!  "filter  capacitance" 
!  "filter  inductance" 

!  "source  resistance" 


! "state  variable  initial  conditions" 

CONSTANT  iLic  =  0.0  ! "Buck  inductor  current" 
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CONSTANT  Voic  =0.0  ! "Buck  output  voltage" 

CONSTANT  Vcfic  =  1100.0  ! "filter  output  voltage" 

CONSTANT  iLfic  =  0.0  ! "filter  inductor  current" 

! "Load  Parameters" 

CONSTANT  Pmax  =  4.5125  1  "100%  load" 

CONSTANT  Rmin  =  45.125  ! "min  load" 

! "start“up" 

CONSTANT  trampup  =  0.1  !  "time  to  ramp  from  0  to  vrefmax" 

CONSTANT  vrefmax  =  950.0  ! "max  ref  output  voltage" 

LOGICAL  rampon 
rampon  =  .TRUE. 

1 "define  the  input  voltage" 

CONSTANT  Vin  =  1100.0 

END  !"of  initial" 

DYNAMIC 

TERMT  (t  .GE.  (tstop~0 . 5*cint ) ) 

DERIVATIVE 

! "Set  the  load  resistance" 

Rload  =  Rmax 

SCHEDULE  rampoff  .XP.  t-trampup 
IF  (rampon)  THEN 

Vref  =  Vrefmax^t /trampup 

ELSE 

Vref  =  Vrefmax 

ENDIF 

!"BUCK  MODELING  EQUATIONS" 

” - determine  conduction  status  of  swl" 

PROCEDURAL  (svif  isw=reframp,  duty,  iL) 

IF  (duty  .LE.  ref ramp)  THEN 

sw  =  .false.  {"switch  open" 

isw  =  0.0  {"source  current  zero" 

ELSE 

sw  =  .true.  {"switch  closed" 
isw  =  iL 

ENDIF 

END  ! "of  procedural" 

”  given  conduction  status  of  sw  and  whether  current  is" 

"  continuous  or  not,  establish  derivative  of  iL" 

PROCEDURAL (piL=sw, contmode, Vcf , Vo, iL, isw) 

if  sw  gated,  assume  can  conduct  pos  current" 

IF  (sw)  THEN 

contmode  =  . true . 

piL  =  (“Rlser*iL+Vcf-Vswdrop“Vo) /Lbuck 
ELSE 

" - if  iL  GT  zero,  diode  conducts" 

IF  (contmode)  THEN 
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piL  =  {-Rlser^iL-Vdidrop“Vo) /Lbuck 
ELSE 

piL  =  0.0  ["discontinuous  mode" 

ENDIF 

ENDIF 

END  l"of  procedural" 

! "link  filter  state  variables  integrated" 

Vcf  =  INTEG( (iLf-isw) /Cf,  Vcfic) 

iLf  =  INTEG( (Vin-Vcf-{Rs*iLf) ) /Lf,  iLfic) 

["establish  the  derivative  of  Buck  output  voltage" 
pVo  =  (iL  -  io)/Cbuck 

["integrate  the  state  variables" 

iLub  =  INTEG(piL,  iLic) 

iL  =  BOUND {-l.Oe-6,  1.0e6,  iLub) 

io  =  Vo/Rload 

Vo  =  INTEG(pVo,  Voic) 

" - form  the  reference  ramp  waveform  (0  to  1)" 

tt  =  mod(t.  Tramp) 
ref ramp  =  tt /Tramp 

" - schedule  discontinuous  mode  when  iL  tries  go  neg" 

SCHEDULE  DCM  .XN.  iL 

[ "Control  Algorithm  for  Buck" 

CONSTANT  hv  =  0.00914 
[CONSTANT  hi  =  -0.003299 
CONSTANT  hi  =  -0.003208 
CONSTANT  hn  =  6.093 

CONSTANT  nic  =  0.0  ' 

CONSTANT  intbnd  =  1000.0 
CONSTANT  minduty  =0.0 
CONSTANT  maxduty  =0.99 
CONSTANT  Kscv  =  1.000 
CONSTANT  Ksci  =  1.000 
va  =  Vref-Vo 
ioave  =  Ksci*io 
iLave  =  Ksci*iL 
ia  =  iLave-ioave 

pn  =  va 

nub  =  INTEG(pn,  nic) 
n  =  BOUND (-intbnd,  intbnd,  nub) 
dss  =  Vref/Vin 

dub  =dss  +  Kscv*hv*va  +  Ksci*hi*ia  +  Kscv*hn*n 
duty  =  BOUND (minduty,  maxduty,  dub) 

END  !  "of  derivative" 

vout  =  Rload*io 

" - establish  transition  to  discontinuous  mode  Buck" 

DISCRETE  DCM 

contmode  =  .false. 


[ "integral  of  volt  err" 

["scale  factor  for  voltage  sensor" 
["scale  factor  for  current  sensor" 
[  "proportional  voltage  term" 

[  "proportional  current  term" 

[  "integral  voltage  term" 
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END 


!"of  dcm 


■  DISCRETE  rampoff 

rampon  =  . FALSE . 

END 

END  !  "of  dynamic" 

END  !  "of  program" 

B.  .CMD  FILE 


! "hardsw.cmd" 

!"Greg  Greseth  -  1  May  1999" 

! "This  file  plots  the  output  for  a  filter-converter  model" 


s  strplt  =  .t. 
s  calplt  =  .f. 


''one  variable  per  x-axis" 


s  devplt  =  1 


s  ppoplt  =  .f. 
s  xinspl  =  6 
s  weditg  =  .f. 
s  wesitg  =  -f. 


!  "1  for  X-windows" 

!  "5  for  postscript" 

I  "true  rotates  plot  90  deg" 

!  "x-axis  plot  units" 

!  "false  suppresses  data  write" 
!  "each  time  SCHEDULE  occurs" 


s  nrwitg  =  .f, 


!  "true  enables  accumulation  of  data" 
1  "after  a  CONTIN" 


s  alcplt  =  .f. 


*false  makes  all  plots  dark" 


prepare  t, Vcf , iLf , iL, io, Vo, duty 

PROCED  startup 
s  ialg=5 
s  maxt=1.0e-6 
s  cint=5.0e-4 
start 

s  pit  =  1 
plot  Vcf,iLf,duty 
s  pit  =  2 
plot  iL,io,Vo 

END 
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APPENDIX  E:  MATLAB  CODE  FOR  IMPEDANCE  CALCULATIONS 


A.  HARD  SWITCHING  MODEL 


%  findxferhard.in 
% 

%  Greg  .Greseth  -  30  March  99 
%  modified  19  April  99 

% 

%  This  program  will  find  the  transfer  function  for  a  linearized  buck 
model 

%  in  order  to  determine  T2  and  T3  for  the  hard  switching  case.  For 
analysis 

%  purposes^,  delta  Vref  is  assumed  to  be  0. 

% 


%  define  the  constants 
Pin  =  200e3; 

Vrefo  =950;  %  =Vout 

Vino  =  1100; 
do  =  0.8636; 

L  =  le-3; 

C  =  2e~3; 

R  =  45.125;  %  =Rcrit 

initdeltaVin  =  5; 

K1  =  -Vrefo/ (Vino) ^2; 
ilo  =  Pin/Vino; 

%  characteristic  equation  coeffiecients 

d2  =  3639.5; 
dl  =  5526354.8; 
d^zero  =  3351129059.4; 

%  calculate  the  gains 

hv  =  ( (L*C) /Vino) * (dl- (1/ {L*C) ) ) ;  %hv  =  0.00914; 

hn  =  (d_zero*L*C) /Vino;  %hn  =  6.093; 

hi  =  -(L/Vino)*(d2-(l/(R*C)));  %hi  =  0.003299; 

Tlnum  =  [hi*C,  hv,  hn] ; 

%  Define  the  A  matrix 

A  =  [ (Vino*hi) /L, (-1/L) - (Vino*hv) /L- (Vino*hi) / (R*L) , (Vino*hn) /L; . . . 
1/C,-1/ (R*C) ,0; . . . 

0,-l,0]; 

%  Define  the  B  matrix 
B  =  [-Vrefo/ (L*Vino) +do/L; 0; 0] ; 
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%  Define  the  D  matrix,  which  is  always  0 
D  =  0; 

%  Define  C  matrix  in  order  to  find  T2 
C_T2  =  [1,0,0]; 

%  Find  the  transfer  function  which  will  define  T2 
[T2niam,  T2den]  =  ss2tf  (A,  B,  C_T2,  D) 

I  Redefine  C  matrix  in  order  to  find  T3 
C_T3  =  [0,1,0]; 

%  Find  the  transfer  function  which  will  define  T3 
[T3num, T3den]  =  ss2tf (A, B, C_T3, D) 

%  Determine  the  Zin  numerator 
Zin_num_hard  =  [T2den,  0] 

%  Determine  the  parts  of  the  Zin  denomenator 

parti  =  [do*T2num,  0] ; 

part2  =  [ilo*Kl*T3den,  0] ; 

part3  =  ilo*conv (Tlnum,  T3num) ; 

part3  =  part3 (1, 2 : 6) ; 

Zin_den_hard  =  partl+part2+part3 

B.  SOFT  SWITCHING  MODEL 

%  findxfersoft .m  . 

% 

%  Greg  Greseth  -  30  March  99 
%  modified  19  April  99 

% 

%  This  program  will  find  the  transfer  function  for  a  linearized  buck 
model 

%  in  order  to  determine  T2  and  T3  for  the  soft  switching  case .  For 
analysis 

%  purposes,  delta  Vref  is  assumed  to  be  0. 

% 

%  define  the  constants 
Pin  =  200e3; 
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%  =Vout 


Vrefo  =  950; 

Vino  =  1100; 
do  =  0.8636; 

L  =  200e-6; 

C  =  400e-6; 

R  =  45.125;  %  =Rcrit 

initdeltaVin  =  5; 

K1  =  -Vrefo/ (Vino)  ■'2; 
ilo  =  Pin/Vino; 


%  characteristic  equation  coefficients 

d2  =  7299; 

dl  =  22194879.2; 

d_zero  =  27000116071.2; 

%  calculate  the  gains 

hv  =  ( {L*C) /Vino) * (dl-(l/ (L*C) ) ) ;  %hv  =  0.00070508; 

hn  =  (d_zero*L*C) /Vino;  %hn  =  1.9636; 

hi  =  -{L/Vino)*(d2-(l/(R*C)));  %hi  =  0.001317; 

Tlnum  =  [hi*C,  hv,  hn] ; 

%  Define  the  A  matrix 

A  =  [ (Vino*hi) /L, (-1/L) - (Vino*hv) /L- (Vino*hi) / (R*L) , {Vino*hn) /L 
1/C,-1/(R*C) ,0; . . . 

0,-l,0]; 

%  Define  the  B  matrix 
B  =  [-Vrefo/ {L*Vino)+do/L;0;0] ; 

%  Define  the  D  matrix,  which  is  always  0 
D  =  0; 


%  Define  C  matrix  in  order  to  find  T2 
C_T2  =  [1,0,0]; 

%  Find  the  transfer  function  which  will  define  T2 
[T2num,T2den]  =  ss2tf (A, B, C_T2, D) 

%  Redefine  C  matrix  in  order  to  find  T3 
C_T3  =  [0,1,0] ; 

%  Find  the  transfer  function  which  will  define  T3 
[T3num, T3den]  =  ss2tf (A,  B, C_T3, D) 

%  Determine  the  Zin  numerator 
Zin  n\im  soft  =  [T2den,  0] 
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%  Determine  the  parts  of  the  Zin  denomenator 

parti  =  [do*T2num,  0] ; 

part2  =  [ilo*Kl*T3den,  0]  ; 

parts  =  ilo*conv (Tlnum,  T3num) ; 

parts  =  parts (1, 2: 6) ; 

Zin_den_soft  =  partl+part2+part3 


C.  NYQUIST  PLOT  FILE 


%  stability. m 
% 

%  Greg  Greseth  -  19  April  99 

% 

%  This  program  will  find  the  Nyquist  plot  of  the  impedance  ratio 
%  between  the  filter  output  impedance  and  the  converter  input 
%  impedance.  The  desired  result  is  to  have  the  magnitude  of  the 
%  ratio  be  less  than  one  over  all  range  of  angular  frequency. 

%  Define  the  functions  being  analyzed: 

%  Zs  is  output  impedance  of  the  filter 
%  Zi  is  input  impedance  of  the  converter 

%  hard  indicates  values  used  for  hard  switching  converter 
%  soft  indicates  values  used  for  soft  switching  converter 

%  IMPORTANT  !!!!!!!!!!! 

%  Both  findxferhard.m  and  findxfersoft .m  must  be  run  with  the  desired 
%  converter  values  must  be  run  first  to  ensure  that  the  input 
%  impedance  vectors  are  in  the  matlab  variable  list 

%  define  the  filter  parameters  and  source  resistance 

Rs  =  0.01;  %  Source  resistance 

Cf  =  6000e-6;  %  Filter  capacitance  -  8265e“6 

Lf  =  0.5e-3;  %  Filter  inductance  -  0.5e-3 

%  calculate  the  filter  output  impedance 

Zsnum  =  [1/Cf,  Rs/ (Cf*Lf ) ] ;  %  numerator 

Zsden  ==  [1,  Rs/Lf,  1/  (Cf*Lf)  ]  ;  %  denomenator 

%  calculate  the  hard  switching  ratio  for  Nyquist  test 

hardnum  =  conv (Zsnum, Zin_den_hard) ; 
hardden  =  conv (Zsden, Zin_num_hard)  ; 

%  calculate  the  soft  switching  ratio  for  Nyquist  test 
softnum  =  conv (Zsnum, Zin_den  soft); 
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softden  =  conv  (Zsden,  Zin_num_soft)  ; 

%  Nyquist  plot  for  hard  switching  case 
figure  (1) 

nyquist (hardnum, hardden) 

%  Nyquist  plot  for  soft  switching  case 
figure  (2) 

nyquist (softnum, softden) 
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APPENDIX  F:  ACSL  MODEL  FOR  LYAPUNOV-BASED  STABILITY 


A.  .CSL  FILE 


PROGRAM  ! "oneload. csl" 

! "Greg  Greseth  -  13  May  1999" 

!"This  file  finds  the  load  current  and  load  voltage  for  a  single  load 
system" 

!"for  the  Lyapunov-based  stability  analysis.  The  non-linear  system 
equations" 

i"for  the  simplified  model  are  utilized.  The  simulation  event  models  a" 
1 "disturbance  (drop)  in  the  source  voltage." 

INITIAL 

CONSTANT  R=0.01 
CONSTANT  Ll=0.0005 
CONSTANT  Cl=. 008265 
CONSTANT  Pl=200000.0 
ALGORITHM  ialg  =  5 
CONSTANT  tstop  =1.0 
NSTEPS  nstp=l 
MAXTERVAL  maxt  =  l.Oe-4 
MINTERVAL  mint  =  l.Oe-7 
CINTERVAL  cint  =  4.0e-3 
CONSTANT  DROP  =55.0 
CONSTANT  VDCSS  =  1100.0 

END  ! "of  initial" 

DYNAMIC 

DERIVATIVE 

VDC=VDCSS-DROP*step (0. 5)  ! "supply  voltage" 

! "initial  conditions" 

CONSTANT  I1IC=190.0 

CONSTANT  V1IC=1000.0 


"resistance  value" 

"inductor  1  value" 

"capacitor  1  value" 

"load  1  power" 

"Runge-Kutta  4th  order" 

"stop  time  for  integration" 

"max  integration  step  size" 
"min  integration  step  size" 
"data  communication  interval" 
"source  voltage  disturbance" 
"source  voltage  steady  state" 


1  "current" 


I1=INTEG {-R/Ll* (II) -Vl/Ll+VDC/Ll, IlIC) 

!  "conditions  for  VI" 

IF (VI. LE. 0.01)  THEN 
pVl=Il/Cl 

ELSE 
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PVl-Il/Cl-Pl/ (C1*V1) 


ENDIF 

! "voltages" 

Vl=INTEG(pVl,VlIC) 

TERMT(T  .GE  .TSTOP,  'Termination  on  time  limit*) 
END  1 "of  derivative” 

END  ! "of  dynamic" 

END  ! "of  program" 


B.  .CMD  FILE 


1 "oneload.cmd” 

!"Greg  Greseth  -  13  May  1999" 

! "This  file  plots  the  load  current  and  load  voltage  for  a  single  load 
system" 

!"for  the  Lyapunov-based  stability  analysis" 

s  strplt  =  .t.  ! "one  variable  per  x-axis" 

s  calplt  =  .f. 

s  devplt  -  1  1"1  for  x-windows" 

! "5  for  postscript" 

prepare  t,  II,  VI 
procedure  plotl 
start 

plot  II  /xtag=*time*  /tag=*Load  1  Current*  VI  /tag=*Load  1 
Voltage  * 

end 
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APPENDIX  G:  ACSL  MODEL  FOR  INTEGRATED  SYSTEM 


A.  .CSL  FILE 


PROGRAM  ! "system. csl" 

!"Greg  Greseth,  LT,  USN  -  13  May  99" 

! "Last  Revised  -  26  May  99 

! "This  ACSL  file  simulates  a  reduced  order  model  for  a 

synchronous" 

! "machine  and  three  phase  rectifier,  with  the  output  through 

a  DC" 

! "link  filter  and  a  buck  chopper  with  a  resistive  load." 

! "Necessary  inputs/known  values  are:" 

! "Output  from  the  model  is:" 


INITIAL  ! "Section  in  which  constants  and  program 

parameters  are  defined" 


ALGORITHM 

ialg  =  5 

! "Runge-Kutta  4th  order 

CONSTANT  tstop  =0.25 

!”stop  time  for 

integration 

ir 

MAXTERVAL 

maxt  ==  S.e-T 

! "max  integration  step 

size" 

MINTERVAL 

mint  =  l.e-6 

! "min  integration  step 

size" 

CINTERVAL 

cint  =  5.e-4 

! "data  communication 

interval” 

CONSTANT  wb  =  377. 

CONSTANT  rkqpri  =  0.0613 
CONSTANT  rkdpri  =0.24 
CONSTANT  rfdpri  =0.00111 
CONSTANT  J  =  1. 

CONSTANT  T1  =  1. 

CONSTANT  B  =  0. 

CONSTANT  P  =  2.  ! "number  of  poles" 

CONSTANT  Xls  =0.08 

CONSTANT  Xmq  =1.00 

CONSTANT  Xlkqpri  =  0.330 

CONSTANT  Xmd  =1.77 

CONSTANT  Xlfdpri  =  0.137 

CONSTANT  Xlkdpri  =  0.334 

CONSTANT  rs  =  0.00515 

CONSTANT  exfdo  =  0.0042 

CONSTANT  alfahat  =0.0 
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CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 

CONSTANT 


PI  =  3.141593 
radl20  =  2.094395 
rad60  -  1.047198 
rad30  -  0.523599 
Vpkstar  =  1.654 
confact  =  400.0 
sterm  =  0.002 


call" 


u  =  0.05  ! "initialize  u  for  use  in  first  findu  function 


{"define  the  subtransient  reactances" 

Xqpripri  =  Xls+ (Xmq*Xlkqpri) / (Xlkqpri+Xmq) 
Xdpripri  = 

Xls+  {Xmd*Xlfdpri*Xlkdpri)  /  (Xmd*Xlfdpri+Xmd*Xlkdpri& 
+X1 kdpr i *X1 f dpr i ) 


{"Buck  Parameters" 
CONSTANT  Lbuck  =  l.Oe-3 
CONSTANT  Riser  =  0.0 
CONSTANT  Vdidrop  =0.0 
CONSTANT  Vswdrop  =0.0 
CONSTANT  Cbuck  =  2 . Oe-3 
CONSTANT  Tramp  =  250.0e-6 


"main  inductance" 

"Lbuck  series  resistance" 
"diode  forw  conduction  drop” 
"switch  forward  drop" 

"main  output  capacitance" 
"switching  period  (4kHz)" 


I  »i - define  and  initialize  the  main  switch  conduction  status" 

LOGICAL  sw 

sw  =  .false.  { "T  sw  conducts,  F  sw  blocks" 

I  - define  and  initialize  variable  which  keeps  track  of" 

{"  whether  the  converter  is  in  contin  or  discontin  mode” 
LOGICAL  contmode 
contmode  =  . true . 

{ "link  filter  parameters" 

CONSTANT  Cf  =  8265. Oe-6  {"filter  capacitance” 

CONSTANT  Lf  =  0.5e“3  {"filter  inductance” 

CONSTANT  Rsrc  =  0.01  {"source  resistance" 


{ "state  variable  initial  conditions" 

CONSTANT  iLic  =  0.0  {"Buck  inductor  current” 

CONSTANT  Voic  =0.0  {"Buck  output  voltage" 

CONSTANT  Vcfic  =  1100.0  {"filter  output  voltage" 
CONSTANT  iLfic  =  0.0  {"filter  inductor  current" 


{"Load  Parameters" 

CONSTANT  Rmax  =  4.5125  {"100%  load" 

CONSTANT  Rmin  =  45.125  { "min  load" 

{ "start-up" 

CONSTANT  trampup 
CONSTANT  vrefmax 
LOGICAL  rampon 
rampon  =  . TRUE . 

{"scale  between  per  unit  values  and  actual  values" 


0.1  {"time  to  ramp  from  0  to  vrefmax* 

950.0  {"max  ref  output  voltage" 
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CONSTANT  currscal  =  202.0 
CONSTANT  volt seal  =  411.0 


COMMON  /  forsub  / 

radSO,  rad60,  radl20,Xqpripri,Xdpripri,  PI,  currscal 


END  ! "of  initial  section” 

DYNAMIC 

TERMT (t.ge. (tstop“0.5*cint) ,  'Termination  on  time  limit*) 
DERIVATIVE 

! "equations  included  for  field  voltage  control  of  the 

generator” 

Vfdpri  =  (exfd^rfdpri) /Xmd 
Vpk  =  sqrt (Vqsr*Vqsr+Vdsr*Vdsr) 

Vdiff  =  Vpkstar  -  Vpk 

exfdub  =  contact *realpl (sterm, Vdiff , exf do) 
exfd  =  BOUND  (-7.0,  7.0,  exfdub) 

! "state  equations  for  psikqpri,  psikdpri,  psifdpri,  wr 

and  thetar” 

psikqpri  ==  INTEG  (-wb*rkqpri*ikqpri,  -0.120) 

psikdpri  =  INTEG {-wb*rkdpri*ikdpri,  1.647) 

psifdpri  =  INTEG  (-wb^rfdpri'^^ifdpri+wb^Vfdpri,  1.777) 

!"wrm  =  INTEG( ( (l./J) * {Te-Tl-B*wrm) ) ,  0.)" 

! "relationship  between  wrm  and  state  wr 

!"wr  =  {P/2.)*wrm” 
wr  =  377. 

thetar  =  INTEG(wr,  0.) 

! "thetar  =  mod (theta, rad60 )  +  radSOd” 

! "electrical  torque  equation" 

! "Te  =  (  (3 . /P) / (4 . *wb) ) * (psidsr*iqsr-psiqsr*idsr) " 

Te  =  psidsr*iqsr-psiqsr*idsr 

! "solve  first  for  ikdpri,  ifdpri  and  ikqpri" 

ikdpri  =  (1. / ( (Xlkdpri+Xmd) * (Xlfdpri+Xmd) - 
Xmd**2 . ) ) * ( (Xlfdpri+Xmd) & 
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*  (psikdpri-Xind*idsr)  -Xitid*  {psifdpri-Xmd*idsr) ) 

ifdpri  =  (1./ { (Xlkdpri+Xmd) * (Xlfdpri+Xmd)-Xmd**2. ) ) * (- 

Xmd* (psikdpriS 

-Xmd*idsr) + (Xlkdpri+Xmd) * (psifdpri-Xmd*idsr) ) 

ikqpri  =  psikqpri/ (Xlkqpri+Xmq) - 
(Xmq*iqsr) / (Xlkqpri+Xmq) 

("equations  to  eliminate  psiqsr  and  psidsr" 

psiqsr  =  Xqpripri*iqsr+psiqpripri 

psidsr  =  Xdpripri*idsr+psidpripri 

("define  the  subtransient  flux  linkages" 

psiqpripri  =  (Xmq/ (Xmq+Xlkqpri) ) *psikqpri 

psidpripri  = 

( (Xmd*psikdpri) /Xlkdpri+ (Xmd*psifdpri) /Xlfdpri) & 

/ (l.+Xmd/Xlfdpri+Xmd/Xlkdpri) 

("define  Eqpripri  and  Edpripri" 

Eqpripri  =  (wr*psidpripri) /wb 

Edpripri  =  - (wr*psiqpripri) /wb 

( "define  phi  to  use  to  determine  alpha" 

phi  =  AT AN2 (Edpripri, Eqpripri) 

( "call  function  to  determine  beta" 

beta  = 

findbeta (psidpripri, psiqpripri, iLf, alfahat) 

( "define  alpha" 
alpha  =  beta  -  phi 

( "desired  reduced  order  voltage-behind-reactance 

equations" 

Vqsr  =  rs*iqsr+ (wr/wb) *Xdpripri*idsr+Eqpripri 

Vdsr  =  rs*idsr- (wr/wb) *Xqpripri*iqsr+Edpripri 

("equation  for  average  DC  voltage  out  of  rectifier" 

Vdc  =  (3 . *sqrt (3 . ) /PI) * (wr/wb) * (sqrt (psiqpripri**2& 
+psidpripri**2) ) * (cos (alpha) )- (3. /PI) * (wr/wb) * ( ( (Xqpripri+Xdpripri ) /2 . ) & 
* (iLf/currscal) + (Xdpripri-Xqpripri) * (iLf /currscal) *sin (2. *beta+rad30) ) 

("need  to  define  K1  for  use  in  subroutines" 
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Kone  =  sqrt(3.)*{“ 

psiqpripri*sin (beta) +psidpripri*cos (beta) ) & 

+  (  (Xdpripri-Xqpripri) *cos (2 . *beta+radl20) -0.5* (Xqpripri+Xdpripri) ) & 
* (iLf /currscal) 


! "call  function  to  determine  u  using  *half-it’  method” 
u  =  findu (Kone,psiqpripri,psidpripri) 

! "now  with  the  new  u,  we  can  determine  conduction 

currents” 


iqsrcond  =  ( (2 . *sqrt (3 . ) ) /PI ) * (iLf /currscal) & 
* (cos (beta+radl20) -cos (beta+u+rad60) ) 

idsrcond  =  ( (2 . *sqrt (3 . ) ) /PI) * (iLf /currscal) & 
* (sin (beta+radl20) -sin (beta+u+rad60) ) 


method” 


!  "calculate  the  commutating  currents  using  Simpson 


used" 


!  "first,  define  the  5  different  values  of  thetar  to  be 


thetarl  =  beta+rad60 
thetar2  =  beta+rad60+u/4 . 
thetar3  =  beta+rad60+u/2 . 
thetar4  =  beta+rad60+ (3. *u) /4 . 
thetarS  =  beta+rad60+u 

! "call  each  q  and  d  subroutine  5  times  to  calculate 

components" 


PROCEDUR?\L  (Q1  = 

thetarl, Kone, psiqpripri,psidpripri,  iLf) 


CALL 

simpsonq (thetarl, Kone, psiqpripri,psidpripri, iLf  ,Q1) 

END  ! "of  procedural" 

PROCEDURAL  (D1  = 

thetarl, Kone, psiqpripri,psidpripri,  iLf) 

CALL 

simpsond (thetarl, Kone, psiqpripri,psidpripri, iLf,  Dl) 

END  !"of  procedural" 

PROCEDURAL  (Q2  = 

thetar2, Kone,psiqpripri,psidpripri,  iLf) 
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CALL 

simpsonq ( thetar2 , Kone, psiqpripri , psidpripri , iLf , Q2 ) 
END  ! ”of  procedural" 

PROCEDURAL  {D2  = 

thetar2, Kone, psiqpripri, psidpripri, iLf ) 

CALL 

simpsond (thetar2, Kone, psiqpripri, psidpripri, iLf ,D2) 
END  ! "of  procedural” 

PROCEDURAL  (Q3  = 

thetarS, Kone, psiqpripri, psidpripri, iLf ) 

CALL 

simpsonq { thetarS , Kone , psiqpripri , psidpripri ; iLf, Q3 ) 
END  ! "of  procedural" 

PROCEDURAL  {D3  = 

thetarS,  Kone, psiqpripri, psidpripri, iLf ) 

CALL 

simpsond (thetarS, Kone, psiqpripri, psidpripri, iLf, D3) 
END  ! "of  procedural" 

PROCEDURAL  {Q4  = 

thetar 4, Kone, psiqpripri, psidpripri, iLf) 

CALL 

simpsonq (thetar 4 , Kone, psiqpripri, psidpripri , iLf, Q4 ) 
END  ! "of  procedural" 

PROCEDURAL  (D4  = 

thet ar4 , Kone , psiqpripri , psidpripri , iLf ) 

CALL 

simpsond ( thetar4 , Kone, psiqpripri , psidpripri, iLf, D4 ) 
END  ! "of  procedural" 

PROCEDURAL  (Q5  = 

thetarS, Kone, psiqpripri, psidpripri, iLf ) 

CALL 

simpsonq (thetarS, Kone, psiqpripri, psidpripri, iLf, QS ) 
END  1 "of  procedural" 

PROCEDURAL  (DS  = 

thet arS , Kone,  psiqpripri , psidpripri , iLf ) 
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CALL 

simpsonci(thetar5,  Kone,psiqpripri,psidpripri,  iLf ,  D5) 

END  ! "of  procedural” 

! "now  add  all  the  components  together” 
iqsrcom  =  (u/ {PI*4 . ) ) * {Ql+4 . *Q2+2 . *Q3+4 . *Q4+Q5) 
idsrcom  =  (u/ (PI*4 . ) ) * (Dl+4 . *D2+2 . *D3+4 . *D4+D5) 

! "add  the  commutating  and  conduction  components  to  get 

average” 

i "iqsr  and  idsr” 

iqsr  =  iqsrcom+iqsrcond 

idsr  =  idsrcom+idsrcond 

!"Set  the  load  resistance” 

Rload  =  Rmax 


SCHEDULE  rampoff  .XP.  t-trampup 
IF  (rampon)  THEN 

Vref  =  Vrefmax*t/trampup 

ELSE 

Vref  =  Vrefmax 

ENDIF 

!"BUCK  MODELING  EQUATIONS” 

” - determine  conduction  status  of  swl” 

PROCEDURAL (sw, isw=ref ramp, duty, iL) 

IF  (duty  .LE.  ref ramp)  THEN 

sw  =  .false.  ! "switch  open” 

isw  =  0.0  ! "source  current  zero” 

ELSE 

sw  =  .true.  ! "switch  closed” 
isw  =  iL 

ENDIF 

END  ! "of  procedural” 

” - given  conduction  status  of  sw  and  whether  current  is” 

”  continuous  or  not,  establish  derivative  of  iL” 

PROCEDURAL (piL=sw, contmode, Vcf , Vo, iL, isw) 

!” - if  sw  gated,  assume  can  conduct  pos  current” 

IF  (sw)  THEN 

contmode  =  .true. 

piL  -  (-Rlser*iL+Vcf"Vswdrop-Vo) /Lbuck 
ELSE 

” - if  iL  GT  zero,  diode  conducts” 

IF  (contmode)  THEN 

piL  =  ( “Riser *iL~Vdidrop-Vo ) /Lbuck 
ELSE 

piL  =  0.0  ! "discontinuous  mode” 

ENDIF 

ENDIF 
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END 


!  "of  procedural” 


! "link  filter  state  variables  integrated" 

Vcf  -  INTEG( (iLf-isw) /Cf,  Vcfic) 

iLf  -  INTEG( (Vdc*voltscal”Vcf-(Rsrc*iLf) ) /Lf,  iLfic) 

! "establish  the  derivative  of  Buck  output  voltage" 
pVo  =  (iL  -  io)/Cbuck 

{"integrate  the  state  variables" 

iLub  =  INTEG(piL,  iLic) 

iL  =  BOUND (-l.Oe-6,  1.0e6,  iLub) 

io  =  Vo/Rload 

Vo  =  INTEG(pVo,  Voic) 

" - form  the  reference  ramp  waveform  (0  to  1)" 

tt  =  mod(t.  Tramp) 
ref ramp  =  tt /Tramp 

” - schedule  discontinuous  mode  when  iL  tries  go  neg" 

SCHEDULE  DCM  .XN.  iL 


! "Control  Algorithm  for  Buck" 
CONSTANT  hv  =  0.00914 
{CONSTANT  hi  -  -0.003299 
CONSTANT  hi  =  -0.003208 
CONSTANT  hn  =  6.093 


CONSTANT  nic  =0.0 
CONSTANT  intbnd  =  1000.0 
CONSTANT  minduty  =0.0 
CONSTANT  maxduty  =  0.99 
CONSTANT  Kscv  =  1.000 
CONSTANT  Ksci  =  1.000 
va  =  Vref-Vo 
ioave  =  Ksci*io 
iLave  =  Ksci^iL 
ia  =  iLave-ioave 


{  "integral  of  volt  err" 


{ "scale  factor  for  voltage  sensor" 
{  "scale  factor  for  current  sensor" 
{ "proportional  voltage  term" 


{ "proportional  current  term" 


pn  =  va 

nub  =  INTEG(pn/  nic)  {"integral  voltage  term" 

'n  =  BOUND {-intbnd,  intbnd,  nub) 
dss  =  Vref/Vcf 

dub  =dss  +  Kscv*hv*va  +  Ksci*hi*ia  +  Kscv*hn*n 
duty  =  BOUND (minduty,  maxduty,  dub) 


END  { "of  derivative" 


vout  =  Rload*io 

” - establish  transition  to  discontinuous  mode  Buck" 

DISCRETE  DCM 

contmode  =  .false. 

END  {"of  dcm" 


DISCRETE  rampoff 

rampon  =  .FALSE. 

END 
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END  ! ”of  dynamic" 

END  ! "of  program" 

SUBROUTINE  simpsonq { thetarsimp, Kone, psiqpripri, psidpripri,  iLf 
+  /iqsrout) 

**★*★*****★★**★**★*★****★*★★*•*•*■*•****  *★***★*****★*******★***■*•★****  *•*•***** 
* 

*  subroutine  to  calculate  components  for  simpson  method  of  determining 

*  the  q  commutating  currents 

* 

******************************************************************************* 
COMMON  /  forsub  /  rad30, rad60, radl20/ Xqpripri, Xdpripri, PI , 

+  currscal 

* 

*  declare  variables 

* 

real  ias,  parti,  part2,  partS,  Kone,  iqsrout,  iLf,  currscal 

* 

*  break  the  equation  for  ias  down  into  manageable  parts 

★ 

partl”Kone-sqrt (3. ) * (psiqpripri*cos (thetarsimp+rad30) +psidpripri 
+  *sin (thetarsimp+rad30) ) 

* 

part2= (Xqpripri-Xdpripri) *cos (2 . *thetarsimp+rad60) +Xqpripri 
+  +Xdpripri 

* 

part3=( (Xqpripri-Xdpripri) *cos (2 . *thetarsimp+radl20) +0. 5 
+  * (Xqpripri+Xdpripri) ) * (iLf /currscal) 

* 

*  combine  the  3  parts  to  determine  ias 

* 

ias= {partl/part2) - (part3/part2) 

* 

*  now  solve  for  iqsr  at  the  inputed  thetar 

* 

iqsrout= ( (2 . *sqrt (3 . ) ) /3 . ) * {ias*cos (thetarsimp+rad30) - 
+  (iLf /currscal) *sin (thetarsimp) ) 

* 

end 

SUBROUTINE  simpsond { thetarsimp, Kone, psiqpripri, psidpripri, iLf 
+  ,idsrout) 

★ 

*  subroutine  to  calculate  components  for  simpson  method  of  determining 

*  the  d  commutating  currents 

* 

COMMON  /  forsub  /  rad30, rad60, radl20, Xqpripri, Xdpripri, PI , 

+  currscal 

* 

*  declare  variables 

* 

real  ias,  parti,  part2,  part3,  Kone,  idsrout,  iLf,  currscal 

★ 

*  break  the  equation  for  ias  down  into  manageable  parts 

* 
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partl=Kone-sqrt  (3 . )  *  {psiqpripri*cos  (thetarsimp+rad30)  +psicipripri 
+  •*sin  {thetarsimp+rad30)  ) 

* 

part 2= (Xqpripri-Xdpripri) *cos (2 . *thetarsimp+rad60) +Xqpripri 
+  tXdpripri 

★ 

part3=( (Xqpripri-Xdpripri) *cos (2 . *thetarsimp+radl20) +0.5 
+  * (Xqpripri+Xdpripri) ) * {iLf /currscal) 

* 

*  combine  the  3  parts  to  determine  ias 

★ 

ias==  (partl/part2)  -  (part3/part2) 

★ 

*  now  solve  for  idsr  at  the  inputed  thetar 

★ 

idsrout= ( (2 . *sqrt (3 . ) ) /3 . ) * (ias^sin (thetarsimp+rad30) + 

+  (iLf/currscal) *cos (thetarsimp) ) 


REAL  FUNCTION  findu (Kone, psiqpripri^psidpripri, iLf ) 


function  subprogram  to  calculate  u  using  the  ’half-it’  method 


COMMON  /  forsub  /  rad30,  radGO,^  radl20,  Xqpripri,  Xdpripri,  PI, 

+  currscal 

* 

*  declare  variables  to  be  used 

* 

logical  condition 

real  utemp,  fcom,  fpricom,  Kone,  iLf,  ul,  u2,  u,  currscal 

* 

ul=0.0 

u2=rad60 

u=0.0 

fcom=Kone-sqrt  (3. )  *  (-psiqpripri*sin (beta+u)  +psidpripri 
+  *cos (beta+u) ) - ( (Xqpripri-Xdpripri) *cos (2 . *beta+2 . *u 
+  -radl20) +0.5* (Xqpripri+Xdpripri) ) * (iLf/currscal) 

10  if  (ABS(ul“u2)  .gt.  0.001)  then 
u  =  0 . 5* (ul+u2) 

fcom=Kone-sqrt (3. ) * (-psiqpripri*sin (beta+u) +psidpripri 
+  *cos (beta+u) ) - ( (Xqpripri-Xdpripri) *cos (2 . *beta+2 . *u 

+  -radl20) +0. 5* (Xqpripri+Xdpripri) ) * (iLf/currscal) 

if  (fcom  .le.  0.00)  then 
ul=u 
else 
u2=u 
endif 
goto  10 
endif 
findu=u 
END 

REAL  FUNCTION  f indbeta (psidpripri, psiqpripri, iLf ,  alfahat) 


*  function  subprogram  to  calculate  beta 
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COMMON  /  forsub  /  rad30, rad60, radl20, Xqpripri, Xdpripri, PI , 

+  currscal 

•k 

*  declare  variables  to  be  used 

* 

real  betatop,  betabot,  betamin,  alfahat,  fbmin,  iLf,  currscal 

★ 

betatop  =  PI 
betabot  =  -PI 
betamin  =  0.0 

* 

*  calculate  the  function  to  determine  if  beta  is  greater  than  betamin 

* 

fbmin  =  SQRT (3 .)* (psiqpripri*cos (betamin) +psidpripri*sin (betamin 
+  +2 . * (iLf /currscal) * (Xqpripri-Xdpripri) 

+  *sin (2 . *betamin-rad60) 

10  if  (ABS (betatop-betabot)  .gt.  0.001)  then 
betamin  =  (betatop+betabot ) /2 . 

fbmin  =  SQRT (3 .)* (psiqpripri*cos (betamin) +psidpripri 
+  *sin (betamin) ) +2 . * (iLf /currscal) * (Xqpripri-Xdpripri) 

+  *sin ( 2 . *betamin-rad60 ) 

if  (fbmin  .le.  0.0)  then 
betabot  =  betamin 
else 

betatop  =  betamin 
endif 
goto  10 
endif 

>*■ 

*  now  determine  beta  based  on  the  above  determination 

* 

findbeta  =  alfahat  +  betamin 
END 


B.  .CMD  FILE 

s  strplt  =  .t. 
s  calplt  =  . f . 

s  devplt  =  1 

s  ppoplt  =  .f. 
s  xinspl  =  6 
s  weditg  =  .f. 
s  wesitg  =  .f. 
s  nrwitg  =  .f. 


!  "one  variable  per  x-axis" 

!  ”1  for  X- windows" 

!  ”5  for  postscript” 

!  "true  rotates  plot  90  deg" 

!  "x-axis  plot  units" 

!  "false  suppresses  data  write" 

!  "each  time  SCHEDULE  occurs" 

!  "true  enables  accumulation  of  data" 


!  "after  a  CONTIN" 


s  alcplt  =  .f.  1  "false  makes  all  plots  dark” 

s  tstop  =  1.0 
s  maxt=5.0e-7 
s  cint=1.0e-3 


prepare 

t, Vdc, Te,psikqpri,psikdpri,psifdpri, iqsr, idsr, u, thetar^beta, alpha, phi& 
Vcf , iLf , duty, Vo, io, exfd 

!  "runs  a  start  up  to  100%,  downpower  to  10%,  and  uppower  to  100%" 

PROCED  runsystem 

s  tstop  =  0.5 
start 

s  cint  =  le-4 
s  devplt  =  5 
s  pit  =  1 

plot  /xlo=0.0  /xhi=0.5  Vdc,Vcf,Vo 
s  pit  =  2 

plot  /xlo=0.0  /xhi=0.5  iLf,io 
s  pit  =  1 

plot  /xlo=0.0  /xhi=0.5  Vdc,Vcf,Vo 
s  pit  =  2 

plot  /xlo=0.0  /xhi=0.5  iLf,io 
s  devplt  =  1 
s  tstop  =  0.52 
contin  . 

s  nrwitg  =  .t. 
s  Rmax  =  45.125 
s  tstop  =  0.57 
contin 

s  Rmax  =  4 . 5125 
s  tstop  =0.62 
contin 

s  devplt  =  5 
s  pit  =  4 

plot  /xlo=0.5  /xhi=0.62  Vdc,Vcf,Vo 
s  pit  =  10 

plot  / xlo=0 . 5  /xhi=0 . 62  iLf , io 
s  pit  =  4 

plot  /xlo=0.5  /xhi=0.62  Vdc,Vcf,Vo 
s  pit  =  10 

plot  /xlo=0.5  /xhi=0.62  iLf,io 
s  devplt  =  1 
s  nrwitg  =  .f. 
s  tstop  =  0.65 
contin 

s  tstop  =  0.67 
s  Cf  =  400e“6 
contin 

s  nrwitg  =  .t. 
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s  Rmax  =  45 . 125 
s  tstop  =  0.12 
contin 

s  Rmax  =  4 . 5125 


s  tstop  =  0.77 
contin 
s  pit  =  13 
s  devplt  =  5 
plot  /xlo=0.65 
s  pit  =  14 
plot  /xlo=0.65 
s  pit  =  13 
plot  /xlo=0.65 
s  pit  =  14 
plot  /xlo=0.65 

END 


/xhi=0.77  Vdc,Vcf,Vo 
/xhi=0.77  iLf,io 
/xhi=0.77  Vdc,Vcf,Vo 
/xhi=0.77  iLf^io 


!"runs  a  start  up  to  10%,  uppower  to  100%  and  downpower  to  10%" 

PROCED  runsystem2 

s  tstop  =0.5 
s  Rmax  =  45.125 
start 

s  cint  =  le-4 
s  devplt  =  5 
s  pit  =  14 

plot  /xlo=0.0  /xhi=0.5  Vdc,Vcf,Vo 
s  pit  =15 

plot  /xlo=0.0  /xhi=0.5  iLf,io 
s  pit  =  14 

plot  /xlo=0.0  /xhi=0.5  Vdc,Vcf,Vo 
s  pit  =  15 

plot  /xlo=0.0  /xhi=0.5  iLf,io 
s  devplt  =  1 
s  tstop  =  0.52 
contin 

s  nrwitg  =  .t. 
s  Rmax  =  4 . 5125 
s  tstop  =  0.57 
contin 

s  Rmax  =  45.125 
s  tstop  =  0.62 
contin 

s  devplt  =  5 
s  pit  =  16 

plot  /xlo=0.5  /xhi=0.62  Vdc,Vcf,Vo 
s  pit  =  17 

plot  /xlo=0.5  /xhi=0.62  iLf,io 
s  pit  =16 

plot  /xlo=0.5  /xhi=0.62  Vdc,Vcf,Vo 
s  pit  =  17 

plot  /xlo=0.5  /xhi=0.62  iLf,io 
s  devplt  =  1 
s  nrwitg  =  . f . 
s  tstop  =0.65 
contin 
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s  tstop  =  0.67 
s  Cf  =  400e~6 
contin 

s  nrwitg  =  .t. 
s  Rmax  =  4 .5125 
s  tstop  =  0.72 
contin 

s  Rmax  =  45.125 
s  tstop  =  0.77 
contin 
s  pit  =  18 
s  devplt  =  5 

plot  /xlo=0.65  /xhi=0.77  Vdc,Vcf,Vo 
s  pit  =  19 

plot  /xlo=0.65  /xhi=0.77  iLf , io 
s  pit  =  18 

plot  /xlo=0.65  /xhi=0.77  Vdc^Vcf^Vo 
s  pit  =19 

plot  /xlo=0.65  /xhi=0.77  iLf,io 


END 
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